I doubt a closed form expression exists. You could solve numerically, using fzero
fun=@(p) p.^2 .*(1-q).*(2.*q.*(1-p.^2)).^(2.*p) -2.*(1-p).*(1-q.*(1-p)).*(p+q.*(1-p));
p=fzero(fun,[0,1]);
However, it is unclear to me whether multiple solutions for p might exist for a given q, even when both are constrained to [0,1]. If so, fzero will find only one of them.
