clear; % clc;

%    Lab 10 RC charge/discharge
%    our state vector will be (x1,y1,vx1,vy1)


C=1.0;  R=10.0;  V=10.0;
tstart=0; tend=10; dt=3.0;

meth=3;
k=0;
Q=0;
for time=tstart:dt:tend
    k=k+1;
    xx1(k)=time;
    yy1(k)=Q;
    true(k)=C*V*(1.0-exp(-time/(R*C) ) );
    if meth == 1   %  Euler
        Q = Q + dt*RHS(Q,R,C,V);
    end
    if meth == 2  %  predictor corrector
        k1 = dt*RHS(Q,R,C,V);      %  Euler predictor
        k2 = dt*RHS(Q+k1,R,C,V);   %  backward Euler corrector
        Q = Q + 0.5*(k1+k2);
    end
    if meth == 3  % RK4
        k1 = dt*RHS(Q,R,C,V);
        k2 = dt*RHS(Q+0.5*k1,R,C,V);
        k3 = dt*RHS(Q+0.5*k2,R,C,V);
        k4 = dt*RHS(Q+k3,R,C,V);
        Q = Q +(k1+2*k2+2*k3+k4)/6.0;
    end
end

e=100.0*(yy1(k)-true(k))/true(k);  
fprintf('percent error for dt=%f and k=%d meth=%d --> %15.10f %%\n',dt,k,meth,e);

%plot(xx1,true,xx1,yy1);

function dydt=RHS(y,R,C,V)
       dydt= V/R - y/(R*C);
end
