# What CAS / solver / math. software to learn & use for PDE systems arising from LQ dynamic optimization?

So I need to cave in and learn some software package ... Until now I have been a paper and pen-mathematician and nearly managed to avoid mathematical software beyond free Wolfram Alpha. To my surprise, I don't find pieces of code snips or instructions that are close enough to my application, so I want to make sure I choose an appropriate tool before spending weeks to learn it. In both the below problem types, everything will have parameters, but not too many in total, so a graphical representation "with a slider" would certainly be helpful. I can use SymPy, Julia seems to have some development momentum, I do have access to some Maple and Matlab versions (but no Mathematica, although I would go for it if that is the way).

What I need this for, is hardly to be considered too advanced mathematics, though there are systems of nonlinear (no worse than quadratic) PDEs. It is dynamic optimization (discrete and/or continuous time, feedback-form), where all optimization is quadratic with linear constraints and could be solved out symbolically. To my surprise, there seems to be hard to find the right words to google for code for dynamic LQ programs; the feedback-form on state is a must for me, open-loop is insufficient. So I need something that can do the following two++ things:

(arising from coupled dynamic programming with quadratic optimization). Type: $\vec u:\mathbf R^n_+\mapsto \mathbf R^n_+$ or $\mapsto \mathbf R^{n+1}_+$ such that each $u_i =$ second-order polynomial in the elements of the Jacobi matrix (requires cubic tensors to write out in (multi)linear algebra language), and with a fairly nice boundary condition $u_i=0$ when $x_i=0$.
It seems that if I introduce $z_i=\sqrt{x_i}$ I can get an expansion: formally write down a MacLaurin series for each function, and match coefficients. Ugly job by hand, I would certainly want something that can do so and check for convergence. Also, I would like concave or convex directions for the solution functions.
I am considering a discrete-time analogue, which for each $t$ will reduce to something like minimizing $\mathbf h\mapsto \mathbf h'\mathbf A_t\mathbf h-2\mathbf b'_t\mathbf h+c_t$ subject to linear (coordinate-wise) inequality constraints $\mathbf M_t\mathbf h\leq \mathbf d_t$ . The $\mathbf A$ are no worse than semidefinite. I thought it would be fairly easy, because in my problems it is not hard to show that I get continuous value functions that piecewise are convex quadratics, but they depend on the parameters and constraints in a not so human-readable manner. Consider even the univariate case: minimizing $\mapsto a_th^2-2b_th)$ subject to, say, $h\leq d_t$ will yield a lot of nested $\max\{$parameter this, parameter that$\}$ expressions, where the coefficients are determined recursively - I was hoping for symbolic expressions for those, hopefully closed-form when remaining horizon (which I can assume finite!) is a fixed number.
I want symbolic expressions for the coefficients (if not in closed-form, then for each time left - I can assume finite horizon in the discrete-time problem), and I need conditions for what $d_t$ are $\leq1$. Also, for when the value function is $C^1$ and for convexity (not only piecewise so).