function fig317

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        
%  where y = (y1, y2, y3 , ..., yn) is an n-vector


% clear all previous variables and plots
clear *
clf

% get(gcf)
set(gcf,'Position', [1067 584 722 900]);
hold on

tmax=160;
N=200;

%  initial values
y0=zeros(4*N-2,1);
y0(2*N-1)=1;

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);  titleX='ode45';
%[t,y] = ode23s(@rhs,[0 tmax],y0);  titleX='ode23s';

nt=size(t);

for it=1:nt(1)
	if t(it)<20
		nt1=it ;
	end;
	if t(it)<40
		nt2=it ;
	end;
	if t(it)<80
		nt3=it ;
	end;
end
nt4=nt(1);

t(nt1)
t(nt2)
t(nt3)
t(nt4)

subaxis(5,1,1,1,'mt',0.005,'mb',0.04,'mr',0.018,'ml',0.08,'p',0.004);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
	jj=jj+1;
	sol(jj)=y(1,j);
end;

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'ko','markersize',2)
grid on
box on
ylabel('y - axis','fontsize',14,'fontweight','bold')
set(gca,'fontsize',14); 
%t0=round(10*t(nt1))/10;
t0=round(t(1));
say=['t = ', num2str(t0)]; 
text(45,0.85,say,'FontSize',14,'FontWeight','bold')
axis([-N N -0.01 1]);

subaxis(5,1,1,2);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
	jj=jj+1;
	sol(jj)=y(nt1,j);
end;

x=linspace(-N+1,N-1,2*N-1);
hold on
box on
plot(x,sol,'--r')
plot(x,sol,'ko','MarkerSize',2)
grid on
ylabel('y - axis','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14); 
%t0=round(10*t(nt1))/10;
t0=round(t(nt1));
say=['t = ', num2str(t0)]; 
text(45,0.15,say,'FontSize',14,'FontWeight','bold')
axis([-N N -0.15 0.2]);

subaxis(5,1,1,3);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
	jj=jj+1;
	sol(jj)=y(nt2,j);
end;

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'ko','MarkerSize',2)
grid on
box on
ylabel('y - axis','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14); 
t0=round(t(nt2));
say=['t = ', num2str(t0)]; 
text(45,0.15,say,'FontSize',14,'FontWeight','bold')
axis([-N N -0.15 0.2]);

subaxis(5,1,1,4);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
	jj=jj+1;
	sol(jj)=y(nt3,j);
end;

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'ko','MarkerSize',2)
grid on
box on
ylabel('y - axis','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14); 
t0=round(t(nt3));
say=['t = ', num2str(t0)]; 
text(45,0.15,say,'FontSize',14,'FontWeight','bold')
axis([-N N -0.15 0.2]);


subaxis(5,1,1,5);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
	jj=jj+1;
	sol(jj)=y(nt4,j);
end;

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'ko','MarkerSize',2)
grid on
box on
t0=round(t(nt4));
say=['t = ', num2str(t0)]; 
text(45,0.15,say,'FontSize',14,'FontWeight','bold')
axis([-N N -0.15 0.2]);


%axis([0 500 -0.2 4]);
% commands to label each axes
xlabel('n - axis','FontSize',14,'FontWeight','bold')
ylabel('y - axis','FontSize',14,'FontWeight','bold')

%set(gca,'xtick',[0 100 200 300 400 500]);

grid on
% command to put legend into plot
%legend(' a = 1',' a = -1','Location','East');

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 

hold off


%  define f(t,y)
function dy=rhs(t,y)
% y(odd) = y_i	y(even) = v_i
N=200;
a=1;
jj=1;
dy=zeros(4*N-2,1);
dy(jj)=y(jj+1);
dy(jj+1)=a*(y(jj+2)-2*y(jj));
for j=2:2*N-2
	jj=jj+2;
	dy(jj)=y(jj+1);
	dy(jj+1)=a*(y(jj+2)-2*y(jj)+y(jj-2));
end;
jj=jj+2;
dy(jj)=y(jj+1);
dy(jj+1)=a*(-2*y(jj)+y(jj-2));
%pause