# Implicit numerical method using NLsolve.jl

I would like to implement `Lobatto IIIA-IIIB` which is an implicit method. To do so, a non-linear equations need to be solved. Then I used `NLsolve.jl` for my problem

``````using NLsolve

function  NuMethod(fV::Function,fT::Function,x₀,y₀,h,N)
x=[] ; push!(x,x₀); y=[] ; push!(y,y₀);
A = [0  0  0 ; 5/24  1/3  -1/24 ; 1/6  2/3  1/6];
B = [1/6  -1/6  0 ; 1/6 1/3  0 ; 1/6  5/6  0];
b = [1/6 , 2/3 , 1/6] ;

g(K,L) = vcat(
K[1] - fV(y₀+h*(B[1,1]*L[1]+B[1,2]*L[2]+B[1,3]*L[3])),
K[2] - fV(y₀+h*(B[2,1]*L[1]+B[2,2]*L[2]+B[2,3]*L[3])),
K[3] - fV(y₀+h*(B[3,1]*L[1]+B[3,2]*L[2]+B[3,3]*L[3])),
L[1] - fT(x₀+h*(A[1,1]*K[1]+A[1,2]*K[2]+A[1,3]*K[3])),
L[2] - fT(x₀+h*(A[2,1]*K[1]+A[2,2]*K[2]+A[2,3]*K[3])),
L[3] - fT(x₀+h*(A[3,1]*K[1]+A[3,2]*K[2]+A[3,3]*K[3])))

g(v) = g(v[1:3],v[4:6]);     #  a function of a single argument

for i in 1:N;
ini_gauss = vcat(fV(y₀),fV(y₀),fV(y₀),fT(x₀),fT(x₀),fT(x₀))
r = nlsolve(g, ini_gauss)
K = r.zero
x₀ = x₀ + h * ( b[1] * K[1] + b[2] * K[2] + b[3] * K[3]) ;
y₀ = y₀ + h * ( b[1] * k[4] + b[2] * K[5] + b[3] * K[6]) ;
push!(x,x₀) ; push!(y,y₀); # stock my solutions
end
return stack(x, dims=1), stack(y, dims=1); # reformulate solutions to vectors
end
``````

I would really like to use this code for any functions `fV` and `fT` single or multidimensional case. However, I always opposite the error that `MethodError: no method matching +(::Vector{Float64}, ::Float64) For element-wise addition, use broadcasting with dot syntax: array .+ scalar`

This is my test probelm

``````function fKeplerT(p) # fKeplerV = vector field to potentioal energy
return Qdot = [p[1] , p[2] , p[3]];
end

function fKeplerV(q)
r = (√(q'*q))^3;
return Pdot = [-q[1]/r , -q[2]/r , -q[3]/r];
end

h=0.01; N=1000;  x₀ = [0.8 , 0.6 , 0 , 0 , 1 , 0.5];

q₀=[x₀[1], x₀[2], x₀[3]];  p₀=[x₀[4], x₀[5], x₀[6]];

q,p = Lobatto_IIIA_IIIB(fKeplerT,fKeplerV,q₀,p₀,h,N);
``````

Note that Labatto tableaus have a lot of structure to them and it’s really not ideal that you solve them like this.

Thanks, but it’s an implicit and this is one way to do that. isn’t?

You get the same error if you just evaluate `g` on `ini_gauss`, so that needs to be corrected.