(* DG.wl -- Find the derivative matrix for Newton's method *)

n=3
Y=Array[y,n]
F=Function[{t,z},{z[[2]]*z[[3]]*Sin[t]-z[[1]]*z[[2]]*z[[3]],
    -z[[1]]*z[[3]]*Sin[t]+(1/20)*z[[1]]*z[[3]],
    z[[1]]^2*z[[2]]-(1/20)*z[[1]]*z[[2]]}]
K1=Array[k,n]
K2=Array[k,n,n+1]
K=Join[K1,K2]

G1=K1-F[t+c[1]*h,Y+h*(a[1,1]*K1+a[1,2]*K2)]
G2=K2-F[t+c[2]*h,Y+h*(a[2,1]*K1+a[2,2]*K2)]
G=Join[G1,G2];

DG=D[G,{K,1}];
DG[[1]][[1]]
DG[[1]][[2]]

mysin=Function[t,"sin("<>ToString[InputForm[t]]<>")"]
Sin[a+b]/.Sin->mysin

For[i=1,i<=6,i++,
    d=InputForm[G[[i]]/.Sin->mysin];
    s=StringForm["G[`1`] = `2`\n",i,d];
    Print[s]
]

fp=OpenWrite["G.out"]
For[i=1,i<=6,i++,
    d=InputForm[G[[i]]/.Sin->mysin];
    s=StringForm["G[`1`] = `2`\n",i,d];
    WriteString[fp,s]
]
Close[fp]

fp=OpenWrite["DG.out"]
For[i=1,i<=6,i++,
    For[j=1,j<=6,j++,
        d=InputForm[DG[[i]][[j]]/.Sin->mysin];
        s=StringForm["DG[`1`, `2`] = `3`\n",i,j,d];
        WriteString[fp,s]
    ]
]
Close[fp]
