I am solving a large system of ODEs with the DifferentialEquation.jl package. Since I am using a number of for loops to set up the RHS to the du, so I would like to see the RHS matrix before solving. Is there a way I can do that?

Lets say I have the following

tspan=[0.0,1.0]
function flux!(du,u,p,t)
for i in 1:n
dummy=0.0
for j in 1:n
dummy=dummy+u[j]
end
du[i]= max(0,dummy)*u[i]+min(0,dummy)*u[i+1]
end
prob = ODEProblem(flux!, u0, tspan)
sol = solve(prob)

So basically I am trying to put a convolution like function on the right hand side, and trying to see if the correct ODE system is generated or not.

One debugging trick here is, because it’s linear, the matrix can be found by calculating the Jacobian. Just do ForwardDiff.jacobian on f and check the matrix.

I want the matrix A, so that I can compute the coefficients by hand and see that they are indeed what they should be. I am really sorry, I am new to Julia and I don’t understand macros.

Thanks for all the help, the issue was when I was writing the flux as a sum on a new line for each component, only the 1st component was being added somehow. So when I changed the coefficient for the 1st component to zero to see how the others affected the solution, the RHS essentially became zero, as the rest of the components were not being added. That’s why I wanted to look at the RHS.

It took me a while to realize that whatever I sum to the flux after putting a newline is not being added to the RHS.