%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.05: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)+3*h/8,y(n)+3*k1/32+9*k2/32);
    k4=h*f(x(n)+12*h/13,y(n)+1932*k1/2197-7200*k2/2197+7296*k3/2197);
    k5=h*f(x(n)+h,y(n)+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
    k6=h*f(x(n)+h/2,y(n)-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/55);
    yd(n+1)=y(n)+(25*k1/216+1408*k3/2565+2197*k4/4104-k5/5);
    y(n+1)=y(n)+(16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55);
    er(n)=y(n+1)-yd(n+1);
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('R-K-Fehlberg Method','fontsize',14,'fontweight','bold');
grid on;
legend('R-K-Fehlberg','Analytic','Location','NorthWest');

text(1/4*xmax,min(y),['Error = ' num2str(er(m-1)) ' at x = ' num2str(x(m))],'fontsize',12,'fontweight','bold');
hold off;