v(nn)=0.9999999999999;

s=zeros(nn+1,1);

s(1)=0;

for ii = 2: nn+1

s(ii)=v(ii-1);

end s

figure;

plot(xm(:,1),xm(:,2),*b-) xlabel(state function xm1)

%Model 1 2.m

function [f,g,xm] = model1 2(um,xinit,par)

nx=par(1);

nu=par(2);

nb=par(3);

ns=par(4); % ns is small interval nn=par(3)*par(4) ; % nn is total interval sm=zeros(nn,1); for ii =1: nb

for jj= 1:ns

sm((ii-1)*ns+jj)=um(ii)/ns;

xm=zeros(nn+1,nx);

xm(1,:)=xinit;

lm=zeros(nn+1,nx);

ma=nx;

t=0;

it=1;

hs=1/nn;

px=model1 3;

xm=nqq(px,nx,nu,nn,xm,ma,t,it,hs,sm,xm,lm,par);

zz=zeros(nn+1,nx);

zz(1,:) = 0;

ma=1;

t=0;

it=1;

hs=1/nn;

x=xm(:,2);

px=model1 4;

jm=nqq(px,nx,nu,nn,zz,ma,t,it,hs,sm,x,lm,par);

%f=jm(nn+1); g(1)=sum (sm)-1; xf=xm(nn+1,1); ed=model1 5;

f=jm(nn+1)+feval(ed,xf,par);

%Model1 3.m compute differential equation

function ff = model1 3(t,it,z,yin,hs,sm,xm,lm,par)

nn=1/hs;

nb=par(3);

ns=par(4);

rr=nn/(3*ns);

p=par(5);

r=par(7);

c=par(9);

T=par(10);

%for i=1:rr,

% us(3*i-2)=0;

% ur(3*i-2)=0;

% us(3*i-1)=0;

% ur(3*i-1)=par(6)/par(7);

% us(3*i)=par(6)/par(7);

% ur(3*i)=0;

%end

for i=1:rr,

us(3*i-2)=par(6)/par(7);

ur(3*i-2)=0;

us(3*i-1)=0;

ur(3*i-1)=par(6)/par(7);

us(3*i)=0;

ur(3*i)=0;

% [A]

% B]

% [C

ts=zeros(nn,1);

if ns > 1

for ii= 1:nb

for jj = 1 :ns

ts(ii)=ts(ii)+sm((ii-1)*ns+jj);

else

for ii= 1: nb

ts(ii)=sm(ii);

pt=nb*ts(floor((it-1)/ns)+1,1);

ff(1)= T*c*([1-ur(1,floor((it-1)/ns)+1)]*r*z(2)-p*z(1))*pt; ff(2)= T*r*z(2)*[ur- (l,flcor ((it-l)/ns)+l)+us(l, floor((it-1)/ns)+1)*(1-z(2)/((1-d) *z(l)))]*pt;

%Model1 4.m

function ff = model1 4(t,it,z,yin,hs,sm,xm,lm,par) nn=1/hs;

fr=nn*mod(t,hs);

xmt=(1-fr)*xm(it)+fr*xm(it+1);

nn=1/hs;

rr=nn/(3*ns);

p=par(5);

r=par(7);

T=par(10);

%for i=1:rr,

% ur(3*i-2)=0;

% ur(3*i-1)=0;

% ur(3*i)=par(6)/par(7);

%end

for i=1:rr,

ur(3*i-2)=0;

ur(3*i-1)=par(6)/par(7); ur(3*i)=0;

umx=zeros(nn,1); umx(2) = sm(1); for ii=3:nn

umx(ii)=umx(ii-1) + sm(ii-1);

qt=nn*s(floor(it),1)*(t-hs*(it-1))+umx(floor(it),1);

ll=exp(-p*qt*T)*[1-ur(1,floor(it))]*r*xmt;

ff=-nn*ll*sm(floor(it),1)*T;

%Model1 5.m

function ff = model1 5(xf,par)

p=par(5);

T=par(10);

ff=-xf*exp(-p*T);

4. Program D: Three Value-Control Model in Chapter 5

%Project3 1.m computation for the integral