Ex1-Runge Kutta-to-h0.001

//solução através do RK, intervalo h=0.001

clc
clear

yanalitico(1)=1;ya(1)=1;yb(1)=1;yc(1)=1;
h=0.001;
t(1)=0;
tfinal=2;
n=tfinal/h;
for k=1:n-1
    t(k+1)=t(k)+h;
//MODO ANALITICO
    yanalitico(k+1)=exp(t(k+1)^4/4-1.5*t(k+1));
//RUNGE KUTTA PRIMEIRA ORDEM
    ya(k+1)=ya(k)+h*(ya(k)*(t(k)^3-1.5));
//RUNGE KUTTA SEGUNDA ORDEM
    yb(k+1)=yb(k)*((2+h*(t(k)^3-1.5))/(2-h*(t(k+1)^3-1.5)));
//RUNGE KUTTA QUARTA ORDEM
    k1=h*yc(k)*((t(k)^3)-1.5);
    k2=h*(yc(k)+(k1/2))*((t(k)+(h/2))^3-1.5);
    k3=h*(yc(k)+(k2/2))*((t(k)+(h/2))^3-1.5);
    k4=h*(yc(k)+k3)*((t(k)+h)^3-1.5);
    yc(k+1)=yc(k)+1/6*(k1+2*k2+2*k3+k4);
end
figure(1)
plot(t,yanalitico,'yo',t,ya,'r',t,yb,'b',t,yc,'g');
legend('Modo Analítico','RK 1','RK 2','RK 4');
title('Exercício 1 - h=0.001')
xlabel('t (s)')
ylabel('y')