Solving six ODEs with 2 set of variables

I have six ODEs as following:

dx1 = - Beta*x1*x2 - (x1*(p1*x1 - p3*x1))/(2*w1);

dx2 = Beta*x1*x2 - (x2*(p2*x2 - p3*x2))/(2*w2);

dx3 = (x1*(p1*x1 - p3*x1))/(2*w1) + (x2*(p2*x2 - p3*x2))/(2*w2);

dp1 = p1*(Beta*x2 + (p1*x1 - p3*x1)/(2*w1)) - Beta*p2*x2 - (p3*(p1*x1 - p3*x1))/(2*w1);

dp2 = Beta*p1*x1 - p2*(Beta*x1 - (p2*x2 - p3*x2)/(2*w2)) - (p3*(p2*x2 - p3*x2))/(2*w2) - 1;

dp3 = 0;

Beta, w1 and w2 are known constants, variable x has 3 dimensions x1, x2, x3, so does p with p1 p2 p3.

And x is my time variant state space representation and p is costate one. So they are actually x(t) and p(t).

I tried to solve it by dsolve, but it says it cannot find explicit solution. So I tried to use ode45 for numerical solution.

But the example I found on Matlab help center doesn't have 2 different variables ODE.

The closest to to mine, is variable x with 2 dimensions x1 and x2. But in my problem, I have another variable p.

Worst of all, it seems I need initial conditions for x1 x2 x3 p1 p2 p3, I do have them for x1(0)=0.9, x2(0)=0.1, x3(0)=0.

But p1 p2, p3 are actually time variant lagrange multiplier, so I don't have initial conditions for them.

I really don't know what to do.

Any tips or help? I will be really appreciated, this is so hard for my small brain.

Thanks in advance