Yes, Z transforms maybe more neat and tidy than my method, but you have to agree, Z transforms are more abstract.

Not necessarily more abstract than working with any other domain transformation, including Fourier transformations.

In fact, if you skip all theoretical derivation of the z transformation, it is a much easier to understand than basically any other transformation. Going from a difference equation to the z domain and back doesn't even require any integrals.

Plus, the z transformation can easily deal with some aspects of signal processing that you'll have a really hard time describing in the "circuit simulation" approach - for example FIR filters (which have some significant advantages over IIR filters and are therefore encountered quite often in signal processing) - they have no "easy" equivalent analog circuit that you could simulate.

Or just try to describe things like a combination (in series or in parallel) of several filters, sampling rate converters, stability analysis, and many other issues, with the "simulated analog circuit" approach - you'll end up with some true nightmares.

The thing is, I have figured out(without using Z transforms) how run 1st or 2nd order (non homogenous) differential equations in my '51,

That doesn't really pertain to the '51 in particular, unless what you mean to say is that you figured out how to do it in C (or '51 assembly) instead of ACSL/Simulink/Fortran or any other language that's actually designed with numeric simulations in mind.

If only it was easy to get this stuff on the market and sell it?

Selling it won't be the problem, but being cheaper (and profitable at the same time) than all of the other people who have already been doing this for quite a while is going to be a problem. Running a PID controller on a '51 is a fairly standard task nowadays.