There is no analytical way to determine the twist for this (or most) load distributions. You must do it either numerically or through guess and check type approaches.
We often think about the lift distribution $l(y)=c_l q c$ by dividing by the dynamic pressure $q=0.5 \rho V^2$ and focusing on the product of the sectional lift coefficient and the chord $c_l c$.
For an elliptical lift distribution, the downwash distribution is uniform. This means that all wing stations 'see' the same angle of attack. Consequently, if we shape the chord elliptically, we can achieve an elliptical distribution with constant $c_l$. Since the downwash (and local aoa) are uniform, this means an un-twisted elliptical chord wing will achieve an elliptical distribution. This is the only easy case.
For any other load distribution, the downwash distribution will be non-uniform and the twist required to get the desired local lift coefficient will be unknown. In fact, it depends on the influence of the rest of the wing -- a full 3d solution.
Fortunately, there are ways. For example, this paper documents a numerical method that can be used to twist any wing to match any target lift distribution. This method has been implemented in a Python script to work with OpenVSP and VSPAERO.