Butcher

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["NumericalDifferentialEquationAnalysis`"];
{p,s}={3,3};
Clear[Eqs,Vars]
Eqs[p_,s_]:=Flatten[
RungeKuttaOrderConditions[p,s,
ButcherRowSum->True,
ButcherSimplify->True, 
RungeKuttaMethod->Explicit]
];
Vars[s_]:=Flatten[{
Table[(Subscript[b, i]),{i,1,s}],
Table[(Subscript[a, i,j]),{i,1,s},{j,1,i-1}],
Table[(Subscript[c, i]),{i,1,s}]
}];
VarsNoC[s_]:=Flatten[{
Table[(Subscript[b, i]),{i,1,s}],
Table[(Subscript[a, i,j]),{i,1,s},{j,1,i-1}]
}];

{p,s}={4,4};
{Timing[fred=FindInstance[Eqs[p,s],Vars[s]];],
Timing[barney=FindInstance[N[Eqs[p,s]],Vars[s]];]}

{p,s}={4,4};
{Timing[fred=Solve[Eqs[p,s],Vars[s]];],
Timing[barney=NSolve[Eqs[p,s],Vars[s]];]}

Here is the FindMinimum computation.

?a

Clear[Objective,Subscript,aa,bb,cc,i,j]
{p,s}={6,9};
Objective=Eqs[p,s]/.{Equal[a_,b_]:>a-b};
Objective=Objective.Objective;
Obj=Objective/.{a->aa,b->bb,c->cc};
vars=Vars[s]/.{a->aa,b->bb,c->cc};
Timing[NMinimize[Obj,vars]]

Clear[Objective,Subscript,aa,bb,cc,i,j]
{p,s}={6,10};
Objective=Eqs[p,s]/.{Equal[a_,b_]:>a-b};
Objective=Objective.Objective;
Obj=Objective/.{a->aa,b->bb,c->cc};
vars=Vars[s]/.{a->aa,b->bb,c->cc};
Timing[NMinimize[Obj,vars]]

Clear[Objective,Subscript,aa,bb,cc,i,j]
{p,s}={7,10};
Objective=Eqs[p,s]/.{Equal[a_,b_]:>a-b};
Objective=Objective.Objective;
Obj=Objective/.{a->aa,b->bb,c->cc};
vars=Vars[s]/.{a->aa,b->bb,c->cc};
Timing[NMinimize[Obj,vars]]

NMinimize

f[y_List]:={{-12,3},{1,-2}}.y
ySol=y/.NDSolve[{
y'[t]==f[y[t]],
y[0]=={1,2}
},
y,
{t, 0 ,TMax}][[1]]
Plot[ySol[t],{t, 0, TMax}]

Clear[Objective,Subscript,aa,bb,cc,i,j]
{p,s}={8,12};
Objective=Eqs[p,s]/.{Equal[a_,b_]:>a-b};
Objective=Objective.Objective;
Obj=Objective/.{a->aa,b->bb,c->cc};
vars=Vars[s]/.{a->aa,b->bb,c->cc};
Timing[NMinimize[Obj,vars]]

No Simplify

Clear[ASol, bSol,cSol]
{p,s}={6,10};
cSol=Union[{0,1},RandomReal[{0,1},s-2]];
cSubs=Table[(Subscript[c, i])->cSol[[i]],{i,Length[cSol]}];
{Length[DeleteCases[Eqs[p,s]/.cSubs,True]], Length[VarsNoC[s]]}

Timing[
SolSub=FindInstance[
Eqs[p,s]/.cSubs,
VarsNoC[s]
];
If[SolSub=={}, Return[Print["No Solution"]],
{ASol,bSol}={Array[a,{s,s}],Array[b,s]}/.SolSub[[1]];
ASol=LowerTriangularize[ASol,-1];
]
];

Length[

ASol=LowerTriangularize[ASol,-1];
MatrixForm[ASol]
bSol
cSol

{0,0.37799561050449726`,0.8464979232622918`,0.952366381908923`}

{p,s}={5,6};
{Timing[fred=Solve[Eqs[p,s],Vars[p,s]];],
Timing[barney=NSolve[Eqs[p,s],Vars[p,s]];]}

ButcherQuadratureConditions[3]

Clear[s]
TableForm[
Table[
RungeKuttaOrderConditions[p],
{p,2,6}],
TableHeadings->Automatic
]

TableForm[
Table[
RungeKuttaOrderConditions[p,ButcherRowSum->True],
{p,2,5}],
TableHeadings->Automatic
]

TableForm[ 
RungeKuttaOrderConditions[6]
]

Table[
TableForm[RungeKuttaOrderCondition,
{p,3,6}]

Clear[p,s]
Eqs[p,s]//TableForm

fred//Length

fred

{Timing[Solve[Eqs,Vars]],
Timing[NSolve[N[Eqs],Vars]]}

ButcherRowConditions[s,p,
RungeKuttaMethod->Explicit]
ButcherColumnConditions[s,p,
RungeKuttaMethod->Explicit]

vars=Flatten[{
Table[(Subscript[b, i]),{i,1,s}],
Table[(Subscript[a, i,j]),{i,1,s},{j,1,i-1}],
Table[(Subscript[c, i]),{i,1,s}]
}]

N[FindInstance[Eqs,vars]]

NMinimize[Eqs.Eqs,vars]

ODE23 Inside Mathematica

Here is an impelementation of the ODE23 inside Mathematica using the plugin framework.  See "tutorial/NDSolvePlugIns" for details

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
ODE23[]["Step"[rhs_, t_, h_, y_, yp_]] :=Module[
{k0,k1,k2,k3,Δy,Error,AbsTol=(10^(-8)),RelTol=(10^(-4)), Safety=0.9},
k0=h yp;
k1=h rhs[t+h/2,y+k0/2];
k2=h rhs[t+3/4h,y+3/4 k1];
k3=h rhs[t+h,y+2/9k0+ 1/3 k1 + 4/9 k2];
Δy=(2k0+3 k1+4  k2)/9;
Error= h*Norm[-(5/72)k0+(1/12)k1+(1/9)k2-(1/8)k3 ];
hNew=Safety*h*((((AbsTol/Error)))^(1/4));
{hNew, Δy, k3/h}];
ODE23[___]["DifferenceOrder"] := 3
ODE23[___]["StepMode"] := Automatic (* Fixed or Automatic *)

xSol=x/. NDSolve[{
x''[t] +Sin[ x[t]] == 0, 
x[0] == 0.99 π, 
x'[0] == 0},
x,
{t,0,TMax=24π},
Method->ODE23 
][[1]];
ts=InterpolatingFunctionCoordinates[xSol][[1]];
ys=InterpolatingFunctionValuesOnGrid[xSol];
tys=Transpose[{ts,ys}];
Plot[xSol[t],{t,0, TMax},
Epilog->Map[Point,tys],
PlotLabel->Length[tys]
]
ListLogPlot[Drop[RotateLeft[ts]-ts,-1],
PlotRange->All]

TR-BDF

TRBDF={(({{0, 0, 0}, {1-(Sqrt[2]), 1-(Sqrt[2]), 0}, {(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), 1-(Sqrt[2])}})),(({{(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), 1-(Sqrt[2])}, {(1/3) (-1+(Sqrt[2])), -(1/3), -(2/3) (-1+(Sqrt[2]))}})),3,"TRBDF"};

TrBdf2[]["Step"[rhs_, t_, h_, y_, yp_]] :=Module[
{n=Length[y],aD = 1-(Sqrt[2]),aO =(1/(2(Sqrt[2]))),Δy,Error,AbsTol=(10^(-8)),RelTol=(10^(-4)), Safety=0.9},
Mat=IdentityMatrix[n]-
k0=h yp;
k1=h rhs[t+h/2,y+k0/2];
k2=h rhs[t+3/4h,y+3/4 k1];
k3=h rhs[t+h,y+2/9k0+ 1/3 k1 + 4/9 k2];
Δy=(2k0+3 k1+4  k2)/9;
Error= h*Norm[-(5/72)k0+(1/12)k1+(1/9)k2-(1/8)k3 ];
hNew=Safety*h*((((AbsTol/Error)))^(1/4));
{hNew, Δy, k3/h}];
ODE23[___]["DifferenceOrder"] := 3
ODE23[___]["StepMode"] := Automatic (* Fixed or Automatic *)

General Implicit RK Scheme Code (revisted)

I set up a general Implicit RK scheme.  It should use a butcher tableau just like the explicit solver and I would like to implement as tests at least the three schemes 237a,b,c on p99.  They have Tableaus (with no error estimates so I put  ∞ in the locations )

RK237A={(({{(1/3), 0}, {1, 0}})),(({{(3/4), (1/4)}, {∞, ∞}})),∞,"Imp_RK237A"};
RK237B={(({{((5-3(Sqrt[2]))/4), ((7-5(Sqrt[2]))/4)}, {((1+(Sqrt[2]))/4), ((3-(Sqrt[2]))/4)}})),(({{((1+(Sqrt[2]))/4), ((3-(Sqrt[2]))/4)}, {∞, ∞}})),∞,"ImplicitRK237B"};
RK237C={(({{(1/4), (1/4)-((Sqrt[3])/6)}, {(1/4)+((Sqrt[3])/6), (1/4)}})),(({{(1/2), (1/2)}, {∞, ∞}})),∞,"ImplicitRK237C"};
TRBDF={(({{0, 0, 0}, {1-(Sqrt[2]), 1-(Sqrt[2]), 0}, {(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), 1-(Sqrt[2])}})),(({{(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), 1-(Sqrt[2])}, {(1/3) (-1+(Sqrt[2])), -(1/3), -(2/3) (-1+(Sqrt[2]))}})),3,"TRBDF"};

Clear[w]
Simplify[(({{0, 0, 0}, {d, d, 0}, {w, w, d}}))//.{γ->2-2(Sqrt[2]),d->γ/2,w->(Sqrt[2])/4}]

Simplify[{w,w,d}-{(1-w)/3,(3w+1)/3,d/3}//.{γ->2-2(Sqrt[2]),d->γ/2,w->(Sqrt[2])/4}]

The first thing to do is work out how to set up the implicit system that needs to be solved.  
For "A" I need to solve
	({{(Subscript[k, 1]), =, f[y+(h/3)(Subscript[k, 1])]}, {(Subscript[k, 2]), =, f[y+h (Subscript[k, 1])]}})
For "B" I need to solve
	({{(Subscript[k, 1]), =, f[y+h ( ((5-3(Sqrt[2]))/4)(Subscript[k, 1])+((7-5(Sqrt[2]))/4) (Subscript[k, 2]))]}, {(Subscript[k, 2]), =, f[y+h ( ((1+(Sqrt[2]))/4)(Subscript[k, 1])+((3-(Sqrt[2]))/4) (Subscript[k, 2]))]}})
For "C" I need to solve
	({{(Subscript[k, 1]), =, f[y+h ( (1/4)(Subscript[k, 1])+((1/4)-((Sqrt[3])/6)) (Subscript[k, 2]))]}, {(Subscript[k, 2]), =, f[y+h (( (1/4)+((Sqrt[3])/6))(Subscript[k, 1])+(1/4) (Subscript[k, 2]))]}})		
"A" is simpler because it only has one implicit equation. "B" and  "C" have two implicit equations.   My iterative fixed point implicit solver needs the A matrix, the step size h, the point y, the function f,  and an initial guess for  ALL the k values.  The solver is going to return the converged values of ALL the k values.   Here is a preliminary version.

ImplicitFixedPoint1[y_,h_,K0_,f_,A_]:= Module[
{s=Length[K0],KNew=K0,i=0, MaxIter=12,yCopies},
yCopies=Table[y,{s}];
Table[
KNew=Map[f,yCopies+h A.KNew], 
{MaxIter}
]
]

Clear[f,k1,k2,h,y,A]
n=4;
f[{y1_,y2_,y3_,y4_}]:= {y1 + (y2^2) Sin[y2 y1],  y4 + Cos[y1 y2],y3 +y2 y1, y1 - y4};
y0=RandomReal[{-1,1},n];
K0=RandomReal[ {-1,1},{2,n}];
A= RK237C[[1]];
h=0.1;
(* For 237A K0=Map[f,{y0+ (1/3)h f[y0],y0+h f[y0]}]; *)
K0=RandomReal[ {-1,1},{2,n}]
Data=ImplicitFixedPoint1[y0,h,K0,f,A];
MatrixForm[Data]

This is seriously inefficient.  What should I do to improve the fixed point iteration?

ImplicitFixedPoint2[y_,h_,K0_,f_,A_]:= Module[
{s=Length[K0],yCopies,KOld=K0,KNew,Tol=1.0 (10^(-7)),i=0, MaxIter=200},
yCopies=Table[y,{s}];
While[
And[
++i<MaxIter,
KNew=Map[f,yCopies+h A.KOld];Norm[KNew-KOld]>Tol],
KOld=KNew
];
Switch[i,
MaxIter,Print["Did not converge with tolerance ", Tol, "!"],
_,Print["Converged after ", i, " iterations!"]
];
KNew
]

Clear[f,k1,k2,h,y,A]
n=4;
f[{y1_,y2_,y3_,y4_}]:= {y1 + (y2^2) Sin[y2 y1],  y4 + Cos[y1 y2],y3 +y2 y1, y1 - y4};
y0=RandomReal[{-1,1},n];
K0=RandomReal[ {-1,1},{2,n}];
A= RK237C[[1]];
h=0.1;
(* For 237A K0=Map[f,{y0+ (1/3)h f[y0],y0+h f[y0]}]; *)
K0=RandomReal[ {-1,1},{2,n}]
Data=ImplicitFixedPoint2[y0,h,K0,f,A];

What would a Newton version look like? The first thing to do is work out what everything is. The equations are 
	({{(Subscript[k, 1]), =, f[y+h (Subscript[a, 1,1])(Subscript[k, 1])+h (Subscript[a, 1,2]) (Subscript[k, 2])]}, {(Subscript[k, 2]), =, f[y+h (Subscript[a, 2,1])(Subscript[k, 1])+h (Subscript[a, 2,2]) (Subscript[k, 2])]}})
and I need to remember that (Subscript[k, i]) are vectors in our vector space (ℝ^n).  Rewriting the equations gives 
	({{(Subscript[g, 1])[(Subscript[k, 1]),(Subscript[k, 2])], :=, (Subscript[k, 1]), -, f[y+h (Subscript[a, 1,1])(Subscript[k, 1])+h (Subscript[a, 1,2]) (Subscript[k, 2])], =, 0}, {(Subscript[g, 2])[(Subscript[k, 1]),(Subscript[k, 2])], :=, (Subscript[k, 2]), -, f[y+h (Subscript[a, 2,1])(Subscript[k, 1])+h (Subscript[a, 2,2]) (Subscript[k, 2])], =, 0}})
and taking the gradient w.r.t (Subscript[k, 1])and (Subscript[k, 2]) gives
	(Subscript[({{(Subscript[∇, 1])(Subscript[g, 1]), =, Id, -, h (Subscript[a, 1,1])∇(Subscript[f, 1])}, {(Subscript[∇, 2])(Subscript[g, 1]), =,  , -, h (Subscript[a, 1,2])∇(Subscript[f, 1])}, {(Subscript[∇, 1])(Subscript[g, 2]), =,  , -, h (Subscript[a, 2,1])∇(Subscript[f, 2])}, {(Subscript[∇, 2])(Subscript[g, 2]), =, Id, -, h (Subscript[a, 2,2])∇f}}), 2])
The structure in the expressions above is worth examining.
*  If (Subscript[a, 11])=(Subscript[a, 22]) then the complicted blocks are similar! Schemes with this characteristic are called Singly Diagonal Implicit
*  As h gets small the matrices become easy to "invert" because the diagonal terms look like the identity and the off-diagonal terms vanish.
* The diagonal matrices are invertible provided h<||(Subscript[a, 1,1])∇(Subscript[f, 1])(Subscript[||, 2]) and h<||(Subscript[a, 2,2])∇(Subscript[f, 2])(Subscript[||, 2]) 
In the fixed point iteration we organized the two vectors into a matrix K={(Subscript[k, 1]),(Subscript[k, 2])}∈(ℝ^(2×n)) .  In the Newton iteration we organize the two vectors into a single long vector 𝓀=(Subscript[k, 1]):(Subscript[k, 2])∈(ℝ^(2n)) and the two vector functions ℊ={(Subscript[g, 1]):(Subscript[g, 2])} into a single long vector because we need to solve linear systems. with this setup our overall gradient becomes
	∇ℊ=(({{Id-h (Subscript[a, 11])∇(Subscript[f, 1]), -h (Subscript[a, 12])∇(Subscript[f, 1])}, {-h (Subscript[a, 21])∇(Subscript[f, 2]), Id-h (Subscript[a, 22])∇(Subscript[f, 2])}}))
and our Newton iteration is
	(Subscript[𝓀, n+1])=(Subscript[𝓀, n])-LinearSolve[∇(Subscript[ℊ, n]),(Subscript[ℊ, n])]
followed by evaluations of ℊ and ∇ℊ.  The complications here are BOOKKEEPING!

ImplicitNewton1[y_,h_,K0_,{f_,df_},A_]:= Module[
{s=Length[K0],n=Length[y],Id,yCopies,KNew=Flatten[K0],MaxIter=12,
g,fder,gDer},
yCopies=Table[y,{s}];
Id=IdentityMatrix[s];
g[k_]:=k- Flatten[Map[f,yCopies+h A.Partition[k,n]]];
fder[k_]:= Map[df, yCopies+h A.Partition[k,n]];
gDer[k_]:=Module[
{dfVals=fder[k]},
IdentityMatrix[s*n]-h ArrayFlatten[
Table[A[[i,j]] dfVals[[j]],{i,s},{j,s}]]
];
Table[
KNew=KNew-LinearSolve[gDer[KNew],g[KNew]],
{MaxIter}
]
]

Clear[f,k1,k2,h,y,A]
n=4;
f[{y1_,y2_,y3_,y4_}]:= {y1 + (y2^2) Sin[y2 y1],  y4 + Cos[y1 y2],y3 +y2 y1, y1 - y4};
df[{y1_,y2_,y3_,y4_}]=D[f[{y1,y2,y3,y4}],{{y1,y2,y3,y4}}];
y0=RandomReal[{-1,1},n];
K0=RandomReal[ {-1,1},{2,n}];
A= RK237C[[1]];
h=0.8;
(* For 237A K0=Map[f,{y0+ (1/3)h f[y0],y0+h f[y0]}]; *)
Fred=ImplicitNewton1[y0,h,K0,{f,df},A]