(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[      8477,        226]
NotebookOptionsPosition[      8233,        213]
NotebookOutlinePosition[      8577,        228]
CellTagsIndexPosition[      8534,        225]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell[BoxData[
 RowBox[{
  RowBox[{"RK23", "[", 
   RowBox[{
   "f_", ",", "u0_", ",", "t0_", ",", "T_", ",", " ", "h0_", ",", " ", 
    "tol_"}], "]"}], ":=", 
  RowBox[{"(", "\[IndentingNewLine]", 
   RowBox[{
    RowBox[{"h", "=", "h0"}], ";", "\[IndentingNewLine]", 
    RowBox[{"t", "=", "t0"}], ";", "\[IndentingNewLine]", 
    RowBox[{"y", "=", 
     RowBox[{"{", "0", "}"}]}], ";", "\[IndentingNewLine]", 
    RowBox[{
     RowBox[{"y", "[", 
      RowBox[{"[", "1", "]"}], "]"}], "=", "u0"}], ";", "\[IndentingNewLine]", 
    RowBox[{"i", "=", "1"}], ";", "\[IndentingNewLine]", 
    RowBox[{"t", "=", 
     RowBox[{"{", "0", "}"}]}], ";", "\[IndentingNewLine]", 
    RowBox[{"While", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"t", "[", 
        RowBox[{"[", "i", "]"}], "]"}], "<", "T"}], ",", 
      "\[IndentingNewLine]", 
      RowBox[{
       RowBox[{"k1", "=", 
        RowBox[{"h", " ", 
         RowBox[{"f", "[", 
          RowBox[{
           RowBox[{"t", "[", 
            RowBox[{"[", "i", "]"}], "]"}], ",", 
           RowBox[{"y", "[", 
            RowBox[{"[", "i", "]"}], "]"}]}], "]"}]}]}], ";", 
       "\[IndentingNewLine]", 
       RowBox[{"k2", "=", 
        RowBox[{"h", " ", 
         RowBox[{"f", "[", 
          RowBox[{
           RowBox[{
            RowBox[{"t", "[", 
             RowBox[{"[", "i", "]"}], "]"}], "+", "h"}], ",", 
           RowBox[{
            RowBox[{"y", "[", 
             RowBox[{"[", "i", "]"}], "]"}], "+", "k1"}]}], "]"}]}]}], ";", 
       "\[IndentingNewLine]", 
       RowBox[{"k3", "=", 
        RowBox[{"h", " ", 
         RowBox[{"f", "[", 
          RowBox[{
           RowBox[{
            RowBox[{"t", "[", 
             RowBox[{"[", "i", "]"}], "]"}], "+", 
            FractionBox["h", "2"]}], ",", 
           RowBox[{
            RowBox[{"y", "[", 
             RowBox[{"[", "i", "]"}], "]"}], "+", 
            RowBox[{
             FractionBox["1", "4"], "k1"}], "+", 
            RowBox[{
             FractionBox["1", "4"], "k2"}]}]}], "]"}]}]}], ";", 
       "\[IndentingNewLine]", 
       RowBox[{"While", "[", 
        RowBox[{
         RowBox[{
          RowBox[{"Abs", "[", 
           RowBox[{
            FractionBox["1", "3"], 
            RowBox[{"(", 
             RowBox[{"k1", "+", "k2", "-", 
              RowBox[{"2", "k3"}]}], ")"}]}], "]"}], ">", "tol"}], ",", 
         "\[IndentingNewLine]", 
         RowBox[{
          RowBox[{"h", "=", 
           RowBox[{"h", 
            SuperscriptBox[
             RowBox[{"(", 
              FractionBox["tol", 
               RowBox[{"Abs", "[", 
                RowBox[{
                 FractionBox["1", "3"], 
                 RowBox[{"(", 
                  RowBox[{"k1", "+", "k2", "-", 
                   RowBox[{"2", "k3"}]}], ")"}]}], "]"}]], ")"}], 
             RowBox[{"1", "/", "3"}]]}]}], ";", "\[IndentingNewLine]", 
          RowBox[{"If", "[", 
           RowBox[{
            RowBox[{"h", "<", 
             SuperscriptBox["10", 
              RowBox[{"-", "5"}]]}], ",", 
            RowBox[{
             RowBox[{
             "Print", "[", 
              "\"\<The step size is too small. Stiff equation suspected\>\"", 
              "]"}], ";", " ", 
             RowBox[{"Return", "[", "]"}]}]}], "]"}], ";", 
          "\[IndentingNewLine]", 
          RowBox[{"k1", "=", 
           RowBox[{"h", " ", 
            RowBox[{"f", "[", 
             RowBox[{
              RowBox[{"t", "[", 
               RowBox[{"[", "i", "]"}], "]"}], ",", 
              RowBox[{"y", "[", 
               RowBox[{"[", "i", "]"}], "]"}]}], "]"}]}]}], ";", 
          "\[IndentingNewLine]", 
          RowBox[{"k2", "=", 
           RowBox[{"h", " ", 
            RowBox[{"f", "[", 
             RowBox[{
              RowBox[{
               RowBox[{"t", "[", 
                RowBox[{"[", "i", "]"}], "]"}], "+", "h"}], ",", 
              RowBox[{
               RowBox[{"y", "[", 
                RowBox[{"[", "i", "]"}], "]"}], "+", "k1"}]}], "]"}]}]}], ";",
           "\[IndentingNewLine]", 
          RowBox[{"k3", "=", 
           RowBox[{"h", " ", 
            RowBox[{"f", "[", 
             RowBox[{
              RowBox[{
               RowBox[{"t", "[", 
                RowBox[{"[", "i", "]"}], "]"}], "+", 
               FractionBox["h", "2"]}], ",", 
              RowBox[{
               RowBox[{"y", "[", 
                RowBox[{"[", "i", "]"}], "]"}], "+", 
               RowBox[{
                FractionBox["1", "4"], "k1"}], "+", 
               RowBox[{
                FractionBox["1", "4"], "k2"}]}]}], "]"}]}]}], ";"}]}], 
        "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", 
       RowBox[{"y", "=", 
        RowBox[{"Append", "[", 
         RowBox[{"y", ",", 
          RowBox[{
           RowBox[{"y", "[", 
            RowBox[{"[", "i", "]"}], "]"}], "+", 
           RowBox[{
            FractionBox["1", "6"], "k1"}], "+", 
           RowBox[{
            FractionBox["1", "6"], "k2"}], "+", 
           RowBox[{
            FractionBox["2", "3"], "k3"}]}]}], "]"}]}], ";", 
       "\[IndentingNewLine]", 
       RowBox[{"t", "=", 
        RowBox[{"Append", "[", 
         RowBox[{"t", ",", 
          RowBox[{
           RowBox[{"t", "[", 
            RowBox[{"[", "i", "]"}], "]"}], "+", "h"}]}], "]"}]}], ";", 
       "\[IndentingNewLine]", 
       RowBox[{"h", "=", 
        RowBox[{"h", 
         SuperscriptBox[
          RowBox[{"(", 
           FractionBox["tol", 
            RowBox[{
             FractionBox["1", "3"], 
             RowBox[{"Abs", "[", 
              RowBox[{"(", 
               RowBox[{"k1", "+", "k2", "-", 
                RowBox[{"2", "k3"}]}], ")"}], "]"}]}]], ")"}], 
          RowBox[{"1", "/", "3"}]]}]}], ";", "\[IndentingNewLine]", 
       RowBox[{"i", "++"}]}]}], "\[IndentingNewLine]", "]"}], ";", 
    "\[IndentingNewLine]", 
    RowBox[{"{", 
     RowBox[{"t", ",", "y"}], "}"}]}], "\[IndentingNewLine]", 
   ")"}]}]], "Input",
 CellChangeTimes->{{3.664871610509589*^9, 3.664871614656519*^9}, {
   3.664871696727519*^9, 3.664871831755891*^9}, {3.6648718619652395`*^9, 
   3.664871954563657*^9}, {3.6648720176121984`*^9, 3.6648720224646025`*^9}, {
   3.6648720694027576`*^9, 3.6648720698631115`*^9}, {3.6648721058325005`*^9, 
   3.664872134340661*^9}, {3.6648721934073696`*^9, 3.664872217575444*^9}, {
   3.6648722598593082`*^9, 3.664872259970392*^9}, {3.6648722916597743`*^9, 
   3.664872305987896*^9}, {3.6648859919958224`*^9, 3.664886012330188*^9}, {
   3.6685006275952377`*^9, 3.6685006998925457`*^9}, {3.6685007378007164`*^9, 
   3.668500927249287*^9}, {3.6685009586664133`*^9, 3.6685009625164323`*^9}, {
   3.6685010818597136`*^9, 3.6685013564686213`*^9}, {3.668501545838112*^9, 
   3.6685015539951315`*^9}, {3.66850278262576*^9, 3.6685028976983414`*^9}, 
   3.668502998349098*^9, {3.6685030341751475`*^9, 3.6685030387584095`*^9}, {
   3.6685031475376315`*^9, 3.668503149301732*^9}, {3.6685031834636865`*^9, 
   3.6685032900097804`*^9}, {3.6685033278939476`*^9, 
   3.6685034179800997`*^9}, {3.6685039380828485`*^9, 
   3.6685039437231708`*^9}, {3.6685048254836044`*^9, 3.668504833079039*^9}, 
   3.668608102483925*^9, {3.668608278807457*^9, 3.6686083806469235`*^9}, {
   3.6686085820739594`*^9, 3.6686085849690523`*^9}, {3.6691016955336423`*^9, 
   3.669101721614134*^9}, {3.669101797545477*^9, 3.6691018039028406`*^9}, 
   3.669101871676717*^9, {3.669101919200435*^9, 3.669101919777468*^9}, {
   3.6691027582564263`*^9, 3.669102787679109*^9}, {3.6691028185438747`*^9, 
   3.669102875848152*^9}, {3.6691047896316147`*^9, 3.6691048618467445`*^9}, 
   3.669104941499301*^9, 3.669105166041144*^9, 3.669108631228341*^9}]
},
WindowSize->{1920, 998},
WindowMargins->{{-8, Automatic}, {Automatic, -8}},
FrontEndVersion->"10.1 for Microsoft Windows (64-bit) (March 23, 2015)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[557, 20, 7672, 191, 739, "Input"]
}
]
*)

(* End of internal cache information *)