Entry 3215

R = 900;
C = 0.001;

h = 0.001;
tmax = 10;
t = (0:h:tmax);

P = zeros(size(t));
P(1) = 0;

PT = zeros(size(t));
PT(1) = 0;

PRK = zeros(size(t));
PRK(1) = 0;

f = 70/60;
Q0 = 20/60;
Q = Q0*sin(2*pi*f.*t);

Q(Q<0)=0;

%% euler
for n = 1:(length(t)-1)
    P(n+1)= P(n) + h*((Q(n)-(P(n)/R))/C);
end

% triprvkovy
for i = 1:length(PT)
    PT(i) = P(i) + 10*Q(i);
end

%% runge kutt

for n = 1:(length(t)-1);
    Qpul =(Q0*sin(2*pi*f*(t(n)+h/2)));
%     Qh = (Q0*sin(2*pi*f*(t(n)+h)));
    Qpul(Qpul<0) = 0;
%     Qh(QQQ<0) = 0;
    
    k1 = (Q(n)-(PRK(n)/R))/C;
    k2 = ((Qpul-((PRK(n)+((h/2)*k1))/R))/C);
    k3 = ((Qpul-((PRK(n)+((h/2)*k2))/R))/C);
    k4 = ((Q(n+1)-((PRK(n)+(h*k3))/R))/C);
    PRK(n+1) = PRK(n)+h/6*(k1+2*k2+2*k3+k4);
end