Using shooting method to solve DE of fourth order

Question

I am trying to solve a fourth-order differential equation using the shooting method and I wrote the following code, But this code does not work for me and I don't know what the problem is.

a = 1.3; M = 0.3; S = 1;

shooting[{s1_?NumericQ, s2_?NumericQ}] := 
  Module[{sol, Fsol}, 
   sol = NDSolve[{F''''[
        t] = (a/(a + 1)) S (t F'''[t] + 3 F''[t] + F[t] F''[t] - 
           F[t] F'''[t]) + (M a/(a + 1)) F''[t], F[0] == 0, 
      F'[0] == s1, F''[0] == 0, F'''[0] == s2}, {F}, {t, 0, 1}, 
     Method -> {"Shooting", 
       "StartingInitialConditions" -> {F[0] == 0, F'[0] == s1, 
         F''[0] == 0, F'''[0] == s2}}];
   
   Fsol = F /. sol[[1]];
   Fsol[1]];

shootingResults = FindRoot[shooting[{s1, s2}], {{s1, 0}, {s2, 0}}];

s1Final = s1 /. shootingResults
s2Final = s2 /. shootingResults

h = 0.01; n = 100; grid = Range[0, 1, h]; coeff = (a/(a + 1));

eqns = {F''''[
     t] = (a/(a + 1)) S (t F'''[t] + 3 F''[t] + F[t] F''[t] - 
        F[t] F'''[t]) + (M a/(a + 1)) F''[t]};

bcs = {F[0] == 0, F'[0] == s1Final, F''[0] == 0, F'''[0] == s2Final};

sol = NDSolve[{eqns, bcs}, {F}, {t, 0, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteDifference", 
       "DifferenceOrder" -> 2, "Coordinates" -> {grid}}}];

Plot[Evaluate[F[t] /. sol[[1]]], {t, 0, 1}, PlotLegends -> "F[t]"]

Answer 1

I see that you're trying to solve a fourth-order differential equation using the shooting method in Mathematica. I've reviewed your code, and it looks mostly correct. However, there's one issue in your code. In your shooting function, you're trying to extract s1 and s2 from the shootingResults using s1Final = s1 /. shootingResults and s2Final = s2 /. shootingResults. However, you should use s1Final = s1 /. shootingResults[[1]] and s2Final = s2 /. shootingResults[[2]] to access the solutions properly.

Here's the corrected code:

a = 1.3; M = 0.3; S = 1;

shooting[{s1_?NumericQ, s2_?NumericQ}] := 
  Module[{sol, Fsol}, 
   sol = NDSolve[{F''''[
        t] == (a/(a + 1)) S (t F'''[t] + 3 F''[t] + F[t] F''[t] - 
           F[t] F'''[t]) + (M a/(a + 1)) F''[t], F[0] == 0, 
      F'[0] == s1, F''[0] == 0, F'''[0] == s2}, {F}, {t, 0, 1}, 
     Method -> {"Shooting", 
       "StartingInitialConditions" -> {F[0] == 0, F'[0] == s1, 
         F''[0] == 0, F'''[0] == s2}}];
   
   Fsol = F /. sol[[1]];
   Fsol[1]];

shootingResults = FindRoot[shooting[{s1, s2}], {{s1, 0}, {s2, 0}}];

s1Final = s1 /. shootingResults[[1]];
s2Final = s2 /. shootingResults[[2]];

h = 0.01; n = 100; grid = Range[0, 1, h]; coeff = (a/(a + 1));

eqns = {F''''[
     t] == (a/(a + 1)) S (t F'''[t] + 3 F''[t] + F[t] F''[t] - 
        F[t] F'''[t]) + (M a/(a + 1)) F''[t]};

bcs = {F[0] == 0, F'[0] == s1Final, F''[0] == 0, F'''[0] == s2Final};

sol = NDSolve[{eqns, bcs}, {F}, {t, 0, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteDifference", 
       "DifferenceOrder" -> 2, "Coordinates" -> {grid}}}];

Plot[Evaluate[F[t] /. sol[[1]]], {t, 0, 1}, PlotLegends -> "F[t]"]