function fig7_15

% requires splineA.m file
%  solve IVP using RK4 then use spline
%  y' = f(t,y)  with   y(0) = y0

M=8;
r=10;
t0=0; y0=0.1;
tmax=1;

% exact solution
nt=100;
tt=linspace(t0,tmax,nt);
a0=(1-y0)/y0;
for it=1:100
    exact(it)=1/(1+a0*exp(-r*tt(it)));
end

tp=linspace(t0,tmax,M+1);
h=tp(2)-tp(1);

y_rk4=rk4('f',tp,y0,h,M+1);
dy0=f(0,y0);
dy1=f(tmax,y_rk4(M+1));
yf = splineA(tp,y_rk4,tt,[1 1],[dy0 dy1]);

clf
% get(gcf)
set(gcf,'Position', [1 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.04,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

hold on
plot(tp,y_rk4,'or','MarkerSize',8,'LineWidth',2)
plot(tt,yf,'b','LineWidth',1.4)
plot(tt,exact,'--k','LineWidth',1.4)
legend({' RK4',' Spline',' Exact'},'Location','SouthEast','FontSize',16,'FontWeight','bold')
xlabel('t-axis')
ylabel('Solution')
axis([0 tmax 0 1.05])
grid on
box on

set(gca,'FontSize',16,'FontWeight','bold')

hold off

%exportgraphics(gcf,'/Users/mark/Desktop/rksp.eps')


% right-hand side of DE
function z=f(t,y)
r=10;
z=r*y*(1-y);

% RK4 method
function ypoints=rk4(f,t,y0,h,n)
y=y0;
ypoints=y0;
for i=2:n
    k1=h*feval(f,t(i-1),y);
    k2=h*feval(f,t(i-1)+0.5*h,y+0.5*k1);
    k3=h*feval(f,t(i-1)+0.5*h,y+0.5*k2);
    k4=h*feval(f,t(i),y+k3);
    yy=y+(k1+2*k2+2*k3+k4)/6;
    ypoints=[ypoints, yy];
    y=yy;
end



