%Differential Equation: dy/dx=y-(1/2)*exp(x/2).*sin(5*x)+5*exp(x/2).*cos(5*x); y(0)=0;
f=@(x,y) y-(1/2)*exp(x/2).*sin(5*x)+5*exp(x/2).*cos(5*x);
f1=@(x,y) exp(x/2).*sin(5*x);   %Analytic Solution
xinit=0; yinit=0;           %y(xinit)=yinit {y(0)=0}

xmax=10;
x=xinit:0.1:xmax;
h=x(2)-x(1);
m=length(x);
y=zeros(1,m);
y(1)=yinit;

for n=1:m-1
    k1=h*f(x(n),y(n));
    k2=h*f(x(n)+h/4,y(n)+k1/4);
    k3=h*f(x(n)+h/4,y(n)+k1/8+k2/8);
    k4=h*f(x(n)+h/2,y(n)-k2/2+k3);
    k5=h*f(x(n)+3*h/4,y(n)+3*k1/16+9*k4/16);
    k6=h*f(x(n)+h,y(n)-3*k1/7+2*k2/7+12*k3/7-12*k4/7+8*k5/7);
    y(n+1)=y(n)+(1/90)*(7*k1+32*k3+12*k4+32*k5+7*k6);
end;

yy=exp(x/2).*sin(5*x);      %Analytic Solution
hold on;
h=plot(x,y,'ro--','LineWidth',1.5);
plot(x,yy,'b-','LineWidth',2.5);
title('5th Order RK Method','fontsize',14,'fontweight','bold');
grid on;
legend('5th Order RK Method','Analytic','Location','NorthWest');
hold off;