(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



Mp[n_,\[Nu]_]:= SparseArray[{
(* Periodic Solution Matrix for CN n is the number of interior pts*)
(* Includes the zeroth point but not the other end *)
Band[{1,1}]->1-2 \[Nu],
Band[{1,2}]->\[Nu],
Band[{2,1}]->\[Nu],
{n+1,1}->\[Nu],
{1,n+1}->\[Nu]
},
{n+1,n+1}]

CNSolver[uIn_List,{\[CapitalDelta]t_, nT_}] :=Module[
{n=Length[uIn]-1, \[CapitalDelta]x,\[Nu],MPlus, MMinus,IntervalLength=2\[Pi],u},
u=uIn;
\[CapitalDelta]x=IntervalLength/(n+1);
\[Nu]=\[CapitalDelta]t/\[CapitalDelta]x^2;
MPlus  =Mp[n,   0.5 \[Nu]];
MMinus=Mp[n,-0.5\[Nu]];
(* Solve MMinus.U^(n+1)= MPlus.U^n repeatedly storing solution in u *) 
Do[
u=LinearSolve[MMinus,MPlus.u],
{nT}];
u
]


uTest1[x_,t_]:=2E^(-3^2t) Sin[3x]+ 4 E^(-2^2t) Cos[2x] ;
nx=96;
\[CapitalDelta]x =2\[Pi]/(nx+1);
u0=Table[uTest1[x,0],{x,0.0,2\[Pi]-\[CapitalDelta]x,\[CapitalDelta]x}] ;
{\[CapitalDelta]t, nt}={0.001,1000};
uEnd=CNSolver[u0,{\[CapitalDelta]t, nt}];
Show[
Plot[uTest1[x, nt*\[CapitalDelta]t],{x,0, 2\[Pi]}],
ListPlot[uEnd,
DataRange->{0, 2\[Pi] - \[CapitalDelta]x}]
]



