How to plot the graph using code?

I have set of partial differential equations, where $F_1=\dfrac{\partial f}{\partial \eta}.$ I am solving them using the NDSolve buil in command of Mathematica. The code for their solution is:

(*Clear all variables*)
ClearAll["Global`*"];
(*Define parameters and constants*)
\[CapitalLambda] = 5;(*Casson fluid parameter*)
Pr = 0.7;(*Prandtl number*)
Sc1 = 160;(*Schmidt number for liquid hydrogen*)
Sc2 = 360;(*Schmidt number for liquid ammonia*)
Le = 10;(*Lewis number*)
Nb = 0.1;(*Brownian diffusion parameter*)
Nt = 0.1;(*Thermophoresis parameter*)
Nr = 0.1;(*Nanoparticles buoyancy ratio parameter*)
Nc1 = 0.1;(*Buoyancy ratio parameters for liquid hydrogen*)
Nc2 = 0.1;(*Buoyancy ratio parameters for liquid ammonia*)
(*\[Chi]=5;Rotation parameter*)
(*\[Lambda]=-2;Mixed convection parameter*)
M = 0.5;

(*Define a range of lambda and chi values*)
\[Lambda]Values = {-2, 0, 5, 
   10};(*Different buoyancy force conditions*)
\[Chi]Values = {0, 2, 5, 10};(*Different rotation strengths*)

(*Define additional terms*)
L1[\[Xi]_] := 1 - Cos[\[Xi]];(*L1 as a function of \[Xi]*)
L2[\[Xi]_] := 1 + Cos[\[Xi]];(*L2 as a function of \[Xi]*)
L3[\[Xi]_] := 2 + Cos[\[Xi]];(*L3 as a function of \[Xi]*)

(*Define coefficients as per the paper*)
B[\[Xi]_] := (1/3)*L3[\[Xi]]*(L2[\[Xi]]^(-1))*
  Tan[\[Xi]/2];(*B as a function of \[Xi]*)
\[Beta]5[\[Xi]_] := (1/2)*L1[\[Xi]]^2*
  L3[\[Xi]];(*beta5 as a function of \[Xi]*)
\[Beta][\[Xi]_] := (2/3)*L3[\[Xi]]*L2[\[Xi]]^2*
  Cos[\[Xi]];(*\[Beta] as a function of \[Xi]*)
s[\[Xi]_] := (8/27)*L3[\[Xi]]*
  L2[\[Xi]]^(-2);(*s as a function of \[Xi]*)
\[Beta]4[\[Xi]_] := 
 L1[\[Xi]]*L2[\[Xi]]^(-1)*
  L3[\[Xi]];(*\[Beta]4 as a function of \[Xi]*)
\[Alpha][\[Xi]_] := \[Chi]*\[Beta][\[Xi]];(*\[Alpha] as a function of \
\[Xi]*)
\[Alpha]1[\[Xi]_] := 
 2*\[Beta][\[Xi]]; (*\[Alpha]1 as a function of \[Xi]*)

(*F1[\[Xi]_,\[Eta]_]:=D[f[\[Xi],\[Eta]],\[Eta]];*)

(*Define system of PDEs*)
eqns =(*Equation of f*){D[
     f[\[Xi], \[Eta]], {\[Eta], 
      3}] + (f[\[Xi], \[Eta]] + 
       2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Xi]]) D[
      f[\[Xi], \[Eta]], {\[Eta], 2}] + \[Beta][\[Xi]] (1 - 
       D[f[\[Xi], \[Eta]], \[Eta]]^2) - 
    2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Eta]] D[
      f[\[Xi], \[Eta]], {\[Xi], 1}, {\[Eta], 
       1}] + \[Lambda] s[\[Xi]] (G[\[Xi], \[Eta]] + 
       Nc1 H[\[Xi], \[Eta]] + Nc2 S[\[Xi], \[Eta]] - 
       Nr P[\[Xi], \[Eta]]) + \[Alpha][\[Xi]] F2[\[Xi], \[Eta]]^2 == 0,
  (*Equation of F2*)
  D[F2[\[Xi], \[Eta]], {\[Eta], 
      2}] + (f[\[Xi], \[Eta]] + 
       2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Xi]]) D[
      F2[\[Xi], \[Eta]], \[Eta]] - \[Alpha]1[\[Xi]] D[
      f[\[Xi], \[Eta]], \[Eta]] F2[\[Xi], \[Eta]] + 
    2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Eta]] D[
      F2[\[Xi], \[Eta]], \[Xi]] == 0,
  (*Equation of G*)
  D[G[\[Xi], \[Eta]], {\[Eta], 2}] + 
    Pr (f[\[Xi], \[Eta]] + 2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Xi]]) D[
      G[\[Xi], \[Eta]], \[Eta]] + 
    Pr (Nb D[P[\[Xi], \[Eta]], \[Eta]] - 
       Nt D[G[\[Xi], \[Eta]], \[Eta]]) D[G[\[Xi], \[Eta]], \[Eta]] - 
    2 B[\[Xi]] Pr D[f[\[Xi], \[Eta]], \[Eta]] D[
      G[\[Xi], \[Eta]], \[Xi]] == 0,
  (*Equation of H*)
  D[H[\[Xi], \[Eta]], {\[Eta], 2}] + 
    Sc1 (f[\[Xi], \[Eta]] + 2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Xi]]) D[
      H[\[Xi], \[Eta]], \[Eta]] - 
    2 B[\[Xi]] Sc1 D[f[\[Xi], \[Eta]], \[Eta]] D[
      H[\[Xi], \[Eta]], \[Xi]] == 0,
  (*Equation of S*)
  D[S[\[Xi], \[Eta]], {\[Eta], 2}] + 
    Sc2 (f[\[Xi], \[Eta]] + 2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Xi]]) D[
      S[\[Xi], \[Eta]], \[Eta]] - 
    2 B[\[Xi]] Sc2 D[f[\[Xi], \[Eta]], \[Eta]] D[
      S[\[Xi], \[Eta]], \[Xi]] == 0,
  (*Equation of F2*)
  D[P[\[Xi], \[Eta]], {\[Eta], 2}] + 
    Le (f[\[Xi], \[Eta]] + 2 B[\[Xi]] D[f[\[Xi], \[Eta]], \[Xi]]) D[
      P[\[Xi], \[Eta]], \[Eta]] - 
    2 B[\[Xi]] Le D[f[\[Xi], \[Eta]], \[Eta]] D[
      P[\[Xi], \[Eta]], \[Xi]] + (Nt/Nb) D[
      G[\[Xi], \[Eta]], {\[Eta], 2}] == 0};

(*Define boundary conditions*)
bcs = {(D[f[\[Xi], \[Eta]], \[Eta]] /. \[Eta] -> 0) == 
    0, (D[f[\[Xi], \[Eta]], \[Eta]] /. \[Eta] -> 5) == 0,
   F2[\[Xi], 0] == 1, F2[\[Xi], 5] == 0,
   G[\[Xi], 0] == 1, G[\[Xi], 5] == 0,
   H[\[Xi], 0] == 1, H[\[Xi], 5] == 0,
   S[\[Xi], 0] == 1, S[\[Xi], 5] == 0,
   P[\[Xi], 0] == 1, P[\[Xi], 5] == 0};

(*Define domain*)
\[Xi]Min = 0; \[Xi]Max = 2;
\[Eta]Min = 0; \[Eta]Max = 5;

(*Solve using NDSolve with MethodOfPDEDiscretization*)
sol = NDSolve[{eqns, bcs}, {f, F2, G, H, S, 
    P}, {\[Xi], \[Xi]Min, \[Xi]Max}, {\[Eta], \[Eta]Min, \[Eta]Max}, 
   Method -> {"PDEDiscretization" -> {"FiniteElement"}}, 
   AccuracyGoal -> 6, PrecisionGoal -> 6, MaxSteps -> 100000];

(*Solve for each parameter combination*)
solutions = 
  Table[NDSolve[{eqns /. {\[Lambda] -> \[Lambda], \[Chi] -> \[Chi]}, 
     bcs}, D[f[\[Xi], \[Eta]], \[Eta]], {\[Xi], \[Xi]Min, \[Xi]Max}, \
{\[Eta], \[Eta]Min, \[Eta]Max}, 
    Method -> {"PDEDiscretization" -> {"FiniteElement"}}, 
    AccuracyGoal -> 6, PrecisionGoal -> 6, 
    MaxSteps -> 
     100000], {\[Lambda], \[Lambda]Values}, {\[Chi], \[Chi]Values}];

How to plot the graph shown using this code.