I am trying to solve the Grad-Shafranov partial differential equation and I am running into the following error

```
using ModelingToolkit
R0 = 2.0
C = 1.0
A = 1.0
#ψ(r,z) = (1/2)*(A*R0^2 +a0*r^2)*z^2 + (1/8)*(C−a0)(r^2−R0^2)^2
@parameters r z
@variables ψ(..)
Dr = Differential(r)
Drr = Differential(r)^2
Dzz = Differential(z)^2
eq = Drr(ψ(r,z)) - Dr(ψ(r,z))/r + Dzz(ψ(r,z)) ~ A*R0^2 + C*r^2
bcs = [ψ(0,z) ~ 0,
ψ(4,z) ~ 0,
ψ(r,-2) ~ 0,
ψ(r,2) ~ 0]
domains = [r ∈ IntervalDomain(0.0,4.0),
z ∈ IntervalDomain(-2.0,2.0)]
pde_system = PDESystem(eq,bcs,domains,[r,z],[ψ])
```

```
MethodError: no method matching ^(::Differential, ::Int64)
Closest candidates are:
^(::Union{AbstractChar, AbstractString}, ::Integer) at strings/basic.jl:718
^(::Union{VectorizationBase.AbstractSIMDVector{W, T}, VectorizationBase.VecUnroll{var"#s2", W, T, V} where {var"#s2", V<:VectorizationBase.AbstractSIMDVector{W, T}}}, ::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}) where {W, T<:Union{Float32, Float64}} at /home/lstagner/.julia/packages/VectorizationBase/qmYqb/src/special/misc.jl:2
^(::Union{VectorizationBase.AbstractSIMDVector{W, T}, VectorizationBase.VecUnroll{var"#s2", W, T, V} where {var"#s2", V<:VectorizationBase.AbstractSIMDVector{W, T}}}, ::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}) where {W, T} at /home/lstagner/.julia/packages/VectorizationBase/qmYqb/src/special/misc.jl:1
...
Stacktrace:
[1] macro expansion
@ ./none:0 [inlined]
[2] literal_pow(f::typeof(^), x::Differential, #unused#::Val{2})
@ Base ./none:0
[3] top-level scope
@ In[4]:11
[4] eval
@ ./boot.jl:360 [inlined]
[5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1094
```

Also, I don’t know how to progress past making the `PDESystem`

All the tutorials were about solving the problem with PINNS but I want to solve it the regular way.