using DifferentialEquations
f = @ode_def Harmonic begin
  dy1 = y2
  dy2 = -y1 - a*y2
end a=>1.5

u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)

functions {
    real[] sho(real t,
    real[] y,
    real[] theta,
    real[] x_r,
    int[] x_i) {
        real dydt[2];
        dydt[1] = y[2];
        dydt[2] = -y[1] - theta[1] * y[2];
    return dydt;
        }
    }
data {
    int<lower=1> T;
    vector[2] y0;
    real ts[T];
    real theta[1];
    }
transformed data {
    real x_r[0];
    int x_i[0];
    }

model {
    }
generated quantities {
    vector[2] y_hat[T];
    matrix[2, 2] A;
    A[1, 1] = 0;
    A[1, 2] = 1;
    A[2, 1] = -1;
    A[2, 2] = -theta[1];
    for (t in 1:T)
    y_hat[t] = matrix_exp((t - 1) * A) * y0;
    // add measurement error
    for (t in 1:T) {
        y_hat[t, 1] = y_hat[t, 1] + normal_rng(0, 0.1);
        y_hat[t, 2] = y_hat[t, 2] + normal_rng(0, 0.1);
        }
    }

println(f.origex)

println(f.params) # prints the parameter of the differential equation