1;

function r=f(x,t);
global n;
nn=n;
La2=2+sin(2*pi*t)+cos(2*pi*t);
La1=1+sin(2*pi*t);
Mu=8100*(1+cos(2*pi*t));
r(1)=-(La1+(nn-2)*La2)*x(1)+Mu*x(2);
for t=2:(nn-1)
r(t)=La2*sum(x(1:(t-2)))+La1*x(t-1)-(La1+(nn-t-1)*La2+Mu)*x(t)+Mu*x(t+1);
endfor
r(nn)=La2*sum(x(1:(nn-2)))+La1*x(nn-1)-Mu*x(nn);
endfunction



global n;
n=200;
t=linspace(0,5,3001);
x0=zeros(n,1); 
x0(1)=1;
y1=lsode("f",x0,t);
ts=linspace(5,6,3001);
y2=lsode("f",y1(end,:),ts);
x0=zeros(n,1); 
x0(n)=1;
y3=lsode("f",x0,t);


sr=0;
sr2=0;
sr3=0;
for k=1:n
sr+=(k-1)*y1(:,k);
sr3+=(k-1)*y3(:,k);
sr2+=(k-1)*y2(:,k);
endfor

filename="Example1"; 


hf = figure ();
hold on;
 plot(t,y1(:,1),'k;p_0|X(0)=0;',t,y3(:,1),'--k;p_0|X(0)=199;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p0.pdf"]);

hf = figure ();
hold on;
 plot(t,y1(:,100),'k;p_99|X(0)=0;',t,y3(:,100),'--k;p_99|X(0)=199;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p99.pdf"]);

hf = figure ();
hold on;
 plot(t,y1(:,200),'k;p_199|X(0)=0;',t,y3(:,200),'--k;p_199|X(0)=199;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p199.pdf"]);


hf = figure ();
hold on;
 plot(ts,y2(:,1),'k;p_0|X(0)=0;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p0S.pdf"]);

hf = figure ();
hold on;
 plot(ts,y2(:,100),'k;p_99|X(0)=0;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p99S.pdf"]);

hf = figure ();
hold on;
 plot(ts,y2(:,200),'k;p_199|X(0)=0;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p199S.pdf"]);

hf = figure ();
hold on;
 plot(t,sr,'k;X(0)=0;',t,sr3,'--k;X(0)=199;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("mean");
grid on;
hold off;
print ([filename,"mean.pdf"]);


hf = figure ();
hold on;
 plot(ts,sr2,'k;X(0)=0;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("mean");
grid on;
hold off;
print ([filename,"meanS.pdf"]);