%
% M-file to do som of the examples in Lewis book.
%

% The frist example is on page 236.

clear all 
A = 1.005; B = 0.01; C=1; P = 3759.9+125; Q = 10; R = 1; N = 100; 

[S10,K10,Kv10] = dricint(A,B,C,P,Q,R,0,N);
Q = 100;
[S100,K100,Kv100] = dricint(A,B,C,P,Q,R,0,N);
Q = 1000;
[S1000,K1000,Kv1000] = dricint(A,B,C,P,Q,R,0,N);

vplot(S10,S100,S1000);
vplot(K10,K100,K1000);
vplot(Kv10,Kv100,Kv1000);

r = vpck([0; 0; 1; 1], [0 50 51 100]); 
k1 = getiv(K10); r = intpol(r,k1);
v10 = dintv(A,B,C,K10,P,10,r);
v100 = dintv(A,B,C,K100,P,100,r);
v1000 = dintv(A,B,C,K1000,P,1000,r);

v10f = dintvf(A,B,C,K10,v10(1),10,r);
v100f = dintvf(A,B,C,K100,v100(1),100,r);
v1000f = dintvf(A,B,C,K1000,v1000(1),1000,r);

vplot(v10,v10f,v100,v100f,v1000,v1000f)

x0 = 1;
[x10,u10,y10] = doptt(A,B,C,K10,Kv10,v10,x0);
[x100,u100,y100] = doptt(A,B,C,K100,Kv100,v100,x0);
[x1000,u1000,y1000] = doptt(A,B,C,K1000,Kv1000,v1000,x0);
vplot(x10,x100,x1000,r);


% look at sub-optimal tracker for Q = 1000;
K = K1000(1); Kv = Kv1000(1); v0 = v1000(1);
[xs,us,ys,vs] = dsoptt(A,B,C,K,Kv,1000,v0,x0,r);
vplot(vs,v1000)
vplot(xs,x1000,r)

Kt = [K1000(1); K1000(1)]; tt = [0 100]; Kt = vpck(Kt,tt);
Kt = intpol(Kt,getiv(r));

v1000t = dintvf(A,B,C,Kt,v1000(1),1000,r);

