1;

function ret=f(x,t);
global n;
La1=10*(2+sin (2*pi*t));
Bs=2+cos (2*pi*t);
ret(1)=-La1*x(1)+Bs*x(n);
for k=2:(n-1)
ret(k)=La1*(x(k-1)-x(k));
endfor
ret(n)=La1*x(n-1)-Bs*x(n);
endfunction

global n;
n=41;

t=linspace(0,14,3001);
x0=zeros(n,1); 
x0(1)=1;
y1=lsode("f",x0,t);
ts=linspace(14,15,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="Example2"; 


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

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

hf = figure ();
hold on;
 plot(t,y1(:,41),'k;p_40|X(0)=0;',t,y3(:,41),'--k;p_40|X(0)=40;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p40.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(:,21),'k;p_20|X(0)=0;')
  legend (hf, "location", "southeast");
 xlabel ("t");
ylabel ("probability");
grid on;
hold off;
print ([filename,"p20S.pdf"]);

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

hf = figure ();
hold on;
 plot(t,sr,'k;X(0)=0;',t,sr3,'--k;X(0)=40;')
  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"]);