1;
function r=f1(x,t)
global n;
la=1+sin(2*pi*t);
mu=3;
r(1)=-la*1.5*x(1)+mu*x(2);
r(2)=la*x(1)-(mu+1.5*la)*x(2)+2*mu*x(3);
for k=3:(n-2)
r(k)=.5*la*x(k-2)+la*x(k-1)-(1.5*la+2*mu)*x(k)+2*mu*x(k+1);
endfor;
r(n-1)=.5*la*x(n-3)+la*x(n-2)-(la+2*mu)*x(n-1)+2*mu*x(n);
r(n)=.5*la*x(n-2)+la*x(n-1)-2*mu*x(n);
endfunction


global n;

h=time();
n=300;
t=linspace(0,100,1000);
tt=linspace(100,101,1000);
x0=zeros(n,1);
x0(1)=1;
y1=lsode("f1",x0,t);
y11=lsode("f1",y1(end,:),tt);
x0=zeros(n,1);
x0(end)=1;
y2=lsode("f1",x0,t);
y22=lsode("f1",y2(end,:),tt);
hh=time()-h;


filename="GBDPEx_"; 
ep=0.015;
epsr=0.015;
hf = figure ();
hold on;
 plot(t,y1(:,1),'k;p_0|X(0)=0;',t,y2(:,1),'--k;p_0|X(0)=299;')
  xlabel ("t");
ylabel ("probability");
grid on;
 axis("tight");
hold off;
print ([filename,"p0.pdf"]);

hf = figure ();
hold on;
plot(tt,y11(:,1),"k",tt,y11(:,1)+ep,"--k",tt,y11(:,1)-ep,"--k");
 xlabel ("t");
ylabel ("probability");
grid on;
 axis("tight");
hold off;
print ([filename,"p0S.pdf"]);

hf = figure ();
hold on;
 plot(t,y1(:,2),'k;p_1|X(0)=0;',t,y2(:,2),'--k;p_1|X(0)=299;')
  xlabel ("t");
ylabel ("probability");
grid on;
 axis("tight");
hold off;
print ([filename,"p1.pdf"]);

hf = figure ();
hold on;
plot(tt,y11(:,2),"k",tt,y11(:,2)+ep,"--k",tt,y11(:,2)-ep,"--k");
 xlabel ("t");
ylabel ("probability");
grid on;
 axis("tight");
hold off;
print ([filename,"p1S.pdf"]);


hf = figure ();
hold on;
 plot(t,y1(:,3),'k;p_2|X(0)=0;',t,y2(:,3),'--k;p_2|X(0)=299;')
  xlabel ("t");
ylabel ("probability");
grid on;
 axis("tight");
hold off;
print ([filename,"p2.pdf"]);

hf = figure ();
hold on;
plot(tt,y11(:,3),"k",tt,y11(:,3)+ep,"--k",tt,y11(:,3)-ep,"--k");
 xlabel ("t");
ylabel ("probability");
grid on;
 axis("tight");
hold off;
print ([filename,"p2S.pdf"]);

hf = figure ();
hold on;
plot(t,y1*(0:(n-1))',"k;X(0)=0;",t,y2*(0:(n-1))',"--k;X(0)=299;")
xlabel ("t");
ylabel ("mean");
grid on;
 axis("tight");
hold off;
print ([filename,"mean.pdf"]);

hf = figure ();
hold on;
%plot(tt,y11*(0:(n-1))',"k;X(0)=0;",tt,y22*(0:(n-1))',"--k;X(0)=20;")
plot(tt,y11*(0:(n-1))',"k",tt,y22*(0:(n-1))'+epsr,"--k",tt,y22*(0:(n-1))'-epsr,"--k");
xlabel ("t");
ylabel ("mean");
grid on;
hold off;
print ([filename,"meanS.pdf"]);

