# Error for solving ODEs by an implicit RK methods when using NLsolve.jl

**URL:** <https://discourse.julialang.org/t/error-for-solving-odes-by-an-implicit-rk-methods-when-using-nlsolve-jl/107764>\
**Category:** General Usage\
**Tags:** question, nlsolve\
**Created:** [December 18, 2023, 10:18am UTC](https://discourse.julialang.org/t/error-for-solving-odes-by-an-implicit-rk-methods-when-using-nlsolve-jl/107764 "2023-12-18T10:18:58Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [December 18, 2023, 10:18am UTC](https://discourse.julialang.org/t/error-for-solving-odes-by-an-implicit-rk-methods-when-using-nlsolve-jl/107764/1 "2023-12-18T10:18:58Z")

</div>

Hi Guys!,

I’m trying to an ODE system using implicit RK method. Defining the following system  
\dot q\_i=p\_i,\quad \dot p\_i = -\frac{q\_i}{||q||^3} and here q,p\in \mathbb R^3. I’m using this code

```julia
using NLsolve

function RK(f::Function,g::Function,t₀,tf,q₀,p₀,N)
    
 h = (tf-t₀)/N;
 t=t₀:h:tf; 

a = [0 0 0 ; 5/24 1/3 -1/24 ; 1/6 2/3 1/6];   
â = [1/6 -1/6 0 ; 1/6 1/3 0 ; 1/6 5/6 0];    

b = [1/6 , 2/3 , 1/6] ;    
    
q=[] ; push!(q,q₀); p=[] ; push!(p,p₀); 

function k_stage!(F,ℓ,q₀,p₀)

F[1] = ℓ[1] .- g(q₀ .+ h * (â[1,1] * ℓ[4] + â[1,2] * ℓ[5] + â[1,3] * ℓ[6]));
F[2] = ℓ[2] .- g(q₀ .+ h * (â[2,1] * ℓ[4] + â[2,2] * ℓ[5] + â[2,3] * ℓ[6]));     
F[3] = ℓ[3] .- g(q₀ .+ h * (â[3,1] * ℓ[4] + â[3,2] * ℓ[5] + â[3,3] * ℓ[6]));
F[4] = ℓ[4] .- f(p₀ .+ h * (a[1,1] * ℓ[1] + a[1,2] * ℓ[2] + a[1,3] * ℓ[3]));
F[5] = ℓ[5] .- f(p₀ .+ h * (a[2,1] * ℓ[1] + a[2,2] * ℓ[2] + a[2,3] * ℓ[3]));        
F[6] = ℓ[6] .- f(p₀ .+ h * (a[3,1] * ℓ[1] + a[3,2] * ℓ[2] + a[3,3] * ℓ[3]));      
        
end

    
for i=1:N;
    
res = nlsolve((F,ℓ) -> k_stage!(F,ℓ,q₀,p₀), rand(6));

k = res.zero
    
q₀ = q₀ + (h/6) * (k[1]+2*k[2]+k[3]);
p₀ = p₀ + (h/6) * (k[4]+2*k[5]+k[6]); 

push!(q,q₀); push!(p,p₀); 
        
end    
  return t, q, p; 
end

```

I opposite an error when applying the above code for the following problem

```julia
function f(p)   
return Qdot = [p[1] , p[2] , p[3]];
end

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

t₀=0; tf=100; 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]]; 

RK(f,g,t₀,tf,q₀,p₀,N)

```

---

<div class="post-metadata">

**Author:** ![JonasWickman](https://avatars.discourse-cdn.com/v4/letter/j/9de0a6/32.png) [@JonasWickman](https://discourse.julialang.org/u/JonasWickman)\
**Post date:** [December 18, 2023, 9:36pm UTC](https://discourse.julialang.org/t/error-for-solving-odes-by-an-implicit-rk-methods-when-using-nlsolve-jl/107764/2 "2023-12-18T21:36:03Z")

</div>

I think the error is because your function `g` returns a 3-vector, which mean that `ℓ[1] .- g(...)` will also be a 3-vector, which you are subsequently trying to store in `F[1]`, which will only store a scalar. I don’t follow what you’re trying to do well enough to suggest solutions though.