v=model1VI1

function v=model1VI1
% MATH LAB RUNGE KUTTA 4 MODEL PACEMAKER VI1 OF ALL THE CARDIAC CELLS WITH
% MATLAB AND SIMULINK
% CARDIAC PARAMETERS
a=-0.3;b=15e-3;do=0.65;d1=-0.0008;R=0.4;E=-90e-3;L=600e-3;C=9600e-6;A=0;F=0;K=0;
% INTIAL CONDITIONS
V(1)=-0.05;I1(1)=0.00096;
% TIME COMPUTATION AND STEP D TIME COMPUTATION,h step
D=2400e-3;N=1000;h=D/N;
for n= 1:1:1000
    A1(n) = (a/C)*(V(n).^3)+(b/C)*(V(n).^2)+(do/C)*V(n)+d1/C -(1/C)*I1(n)+(A*sin(2*(3.14)*(n*h)*F+K))/C;
    B1(n) = (1/L)*V(n)-(R/L)*I1(n)+E/L;
    A2(n) = (a/C)*(V(n)+(h/2)*A1(n)).^3+(b/C)*(V(n)+(h/2)*A1(n)).^2+(do/C)*(V(n)+(h/2)*A1(n))+d1/C-(1/C)*(I1(n)+(h/2)*B1(n))+(A*sin((2*3.14)*F*(n*h)+K))/C;
    B2(n) = (1/L)*(V(n)+(h/2)*A1(n))-(R/L)*(I1(n)+(h/2)*B1(n))+E/L;
    A3(n) = (a/C)*(V(n)+(h/2)*A2(n)).^3 +(b/C)*(V(n)+(h/2)*A2(n)).^2 +(do/C)*(V(n)+(h/2)*A2(n)) +d1/C -(1/C)*(I1(n)+(h/2)*B2(n))+(A*sin((2*3.14)*F*(n*h)+K))/C;
    B3(n) = (1/L)*(V(n)+(h/2)*A2(n))-(R/L)*(I1(n)+(h/2)*B2(n))+E/L;
    A4(n) = (a/C)*(V(n)+h*A3(n)).^3+(b/C)*(V(n)+h*A3(n)).^2+(do/C)*(V(n)+h*A3(n))+d1/C-(1/C)*(I1(n)+h*B3(n))+(A*sin((2*3.14)*F*(n*h)+K))/C;
    B4(n) = (1/L)*(V(n)+h*A3(n))-(R/L)*(I1(n)+h*B3(n))+E/L;
    V(n+1) = V(n)+(h/6)*(A1(n)+A2(n)+2*A3(n)+A4(n));
    I1(n+1) = I1(n)+(h/6)*(B1(n)+B2(n)+2*B3(n)+B4(n));
end
plot(V)
title(['a= ', num2str(a),'  b= ', num2str(b),'    do= ', num2str(do),'    d1= ', num2str(d1),'    E= ', num2str(E),'   L= ', num2str(L),'   C= ', num2str(C),'   D= ', num2str(D),'   N=' , num2str(N),'   R=' , num2str(R)])