Dear Julia Users,

I'm trying to speed up some code that employs passing functions as 
arguments. One part of the code solves an ODE; if I use CVODE from 
Sundials, and rewrite the function to accept callable types, I get a 
ReadOnlyMemoryError - as I don't know the Sundials API, can someone help me 
with where I'm going wrong? Code below.

Best
Simon

**

using Sundials

function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64}; 
reltol::Float64=1e-4, abstol::Float64=1e-6)
    neq = length(y0)
    mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
    flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32, 
(Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})), 
t[1], Sundials.nvector(y0))
    flag = Sundials.CVodeSetUserData(mem, f)
    flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
    flag = Sundials.CVDense(mem, neq)
    yres = zeros(length(t), length(y0))
    yres[1,:] = y0
    y = copy(y0)
    tout = [0.0]
    for k in 2:length(t)
        flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
        yres[k,:] = y
    end
    Sundials.CVodeFree([mem])
    return yres
end

function f(t, y, ydot)
    ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
    ydot[3] = 0.
    ydot[2] = 0.
end

immutable J; end
call(::Type{J},t, y, ydot) = f(t, y, ydot)

t = [0.1:0.0001:1]
res = cvode(J, [-60.0, 0.0, 0.0], t);

Reply via email to