#COPYRIGHT NOTICE:  Copyright (c) 1997-1998 by Sasha Cyganowski.
#Permission is  
#granted to anyone to use, modify, or redistribute this software freely, 
#subject to the following restrictions:  1. The author accepts no
#responsibility for any consequences of this software and makes no guarantee
#that this software is free of defects.  2.  The origin of this software must  
#not be misrepresented, either by explicit claim or by omission. 
#3.  This notice and the copyright must be included with all copies or
#modified versions of this software.  4.  This software may not be included or
#redistributed as part of any package to be sold for profit unlesss the author
#has given explicit written permission to do so.    
# Sasha Cyganowski, School of Computing and Mathematics, Deakin University,
#Waurn Ponds, VIC, Australia.  Internet: sash@deakin.edu.au
# Telephone:  +61 3 5227 1382.

supp1:=proc(a,b,i)
local soln,Y,Delta,n,j,u;
soln[i]:=Y.i[n]+a[i]*Delta[n]+sum('op(j,op(i,b))*Delta*W.j[n]','j'=1..nops(op(1,b)));
for u from 1 to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od;
end:\


supp2a:=proc(a,b,i,j)
local soln,u,Y,n,Delta;
soln[i]:=Y.i[n]+a[i]*Delta[n]+op(j,op(i,b))*sqrt(Delta[n]);

for u from 1 to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od;
end:

supp2b:=proc(a,b,i,j)
local soln,u,Y,n,Delta;
soln[i]:=Y.i[n]+a[i]*Delta[n]-op(j,op(i,b))*sqrt(Delta[n]);

for u from 1 to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od;
end:

supp3a:=proc(b,j,i)
local soln,u,Y,n,Delta;
soln[i]:=Y.i[n]+op(j,op(i,b))*sqrt(Delta[n]);

for u from 1 to nops(b) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od;
end:

supp3b:=proc(b,j,i)

local soln,u,Y,n,Delta;
soln[i]:=Y.i[n]-op(j,op(i,b))*sqrt(Delta[n]);

for u from 1 to nops(b) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od;
end:
supp4a:=proc(a,b,i,j)
local tempb,k;
tempb:=op(j,op(i,b));
for k from 1 to nops(a) do
tempb:=subs(x[k]=supp2a(a,b,k,j),tempb);
od;
end:

supp4b:=proc(a,b,i,j)
local tempc,k;
tempc:=op(j,op(i,b));
for k from 1 to nops(a) do
tempc:=subs(x[k]=supp2b(a,b,k,j),tempc);
od;
end:

supp5a:=proc(b,r,i,j)
local tempd,k;
tempd:=op(j,op(i,b));
for k from 1 to nops(b) do
tempd:=subs(x[k]=supp3a(b,r,k),tempd);
od;
end:

supp5b:=proc(b,r,i,j)
local tempe,k;
tempe:=op(j,op(i,b));
for k from 1 to nops(b) do
tempe:=subs(x[k]=supp3b(b,r,k),tempe);
od;
end:

expwk2:=proc(a,b)
# Derivative free scheme
local tempa,tempb,tempc,tempd,tempe,k,Y,Delta,W,j,i,n,V,soln,r,j1,r1,u,DeltaW;


for i from 1 to nops(a) do
 tempa[i]:=op(i,a):
for k from 1 to nops(a) do
tempa[i]:=subs(x[k]=supp1(a,b,k),tempa[i]);


od;
soln[i]:=Y.i[n+1]=Y.i[n]+1/2*(tempa[i]+op(i,a))*Delta[n]+1/4*sum('((supp4a(a,b,i,j)+supp4b(a,b,i,j)+2*op(j,op(i,b)))*Delta*W.j[n]+(sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*Delta*W.j[n]','r'=1..j-1)+sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*Delta*W.j[n]','r'=j+1..nops(op(1,b)))))','j'=1..nops(op(1,b)))+1/4*sum('((supp4a(a,b,i,j1)-supp4b(a,b,i,j1))*((Delta*W.j1[n])^2-Delta[n])+(sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=1..j1-1)+sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=j1+1..nops(op(1,b)))))*(Delta[n])^(-1/2)','j1'=1..nops(op(1,b)));
for u from 1 to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);

od;
print(soln[i]);
od;
end:

impwk2:=proc(a,b)
# Derivative free scheme
local tempa,k,Y,Delta,W,j,i,n,V,soln,r,j1,r1,u;

for i from 1 to nops(a) do

tempa[i]:=op(i,a):
for k from 1 to nops(a) do
tempa[i]:=subs(x[k]=supp1(a,b,k),tempa[i]);
od;
soln[i]:=Y.i[n+1]=Y.i[n]+1/2*(tempa[i]+op(i,a))*Delta[n]+1/4*sum('(supp4a(a,b,i,j)+supp4b(a,b,i,j)+2*op(j,op(i,b))+(sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*(Delta[n])^(-1/2)','r'=1..j-1)+sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*(Delta[n])^(-1/2)','r'=j+1..nops(op(1,b)))))*Delta*W.j[n]','j'=1..nops(op(1,b)))+1/4*sum('((supp4a(a,b,i,j1)-supp4b(a,b,i,j1))*((Delta*W.j1[n])^2-Delta[n])+(sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=1..j1-1)+sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=j1+1..nops(op(1,b)))))*(Delta[n])^(-1/2)','j1'=1..nops(op(1,b)));
for u from 1 to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);

od;
print(soln[i]);
od;
end: