#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.

LJ:=proc(X,b,j)
local k;
sum('op(j,op(k,b))*diff(X,x[k])','k'=1..nops(b));
end:
L0:=proc(X,a,b)
local part1,part2,part3,k,j,l;
part1:=diff(X,t);
part2:=sum('a[k]*diff(X,x[k])','k'=1..nops(a));
part3:=1/2*sum('sum('sum('op(j,op(k,b))*op(j,op(l,b))*diff(X,x[k],x[l])','j'=1..nops(op(1,b)))','k'=1..nops(a))','l'=1..nops(a));
part1+part2+part3;
end:
WeakImp2:=proc(a,b)
global soln;
local l,u,i,n,Y,Delta,W,j1,j2,V,p,k,temp;
for i from 1 to nops(a) do
temp[i]:=op(i,a):
for k from 1 to nops(a) do
temp[i]:=subs(x[k]=Y.k[n+1],temp[i]);
od;


soln[i]:=Y.i[n+1]=Y.i[n]+1/2*(temp[i]+op(i,a))*Delta[n]+sum('op(j,op(i,b))*Delta*W.j[n]','j'=1..nops(op(1,b)))+1/2*sum('L0(op(l,op(i,b)),a,b)*Delta*W.l[n]*Delta[n]','l'=1..nops(op(1,b)))+1/2*sum('sum('LJ(op(j2,op(i,b)),b,j1)*(Delta*W.j1[n]*Delta*W.j2[n]+V[(j1,j2)])','j1'=1..nops(op(1,b)))','j2'=1..nops(op(1,b)));
for u from 1 to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od;
soln[i]:=Y.i[n+1]=solve(soln[i],Y.i[n+1]);


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