Maple can produce a graphic which is quite closeto the one you display.
But this takes time and it is really for the sake of art (IMO):

Customized_graphic.mw

rt
pattern := "%15s":

for i from 1 to 5 do
  printf(cat(pattern$i, "\n"), "Hello"$i)
end do
          Hello
          Hello          Hello
          Hello          Hello          Hello
          Hello          Hello          Hello          Hello
          Hello          Hello          Hello          Hello          Hello

use the big green arrow in the menubar to upload your worksheet.

@vv 

Here are two files:

• To_vv_1.mw
  I compare the solutions pdsolve returns on your example and a modified one to the solutions I get with the simple approach (SA) I predented in my first reply to @emendes (the one which requires compatibility conditions).
  SA does not provide a solution for your toy problem  (a <> -d as you noticed).
  For the second toy problem the pdsolve and SA give the same fundamental solution U(x, y).
  I name it  fundamental in the sense that any function of U(x, y) is also a solution.
• To_vv_2.mw
  I compared the results from pdsolve to those from the Method Of Characteristics (MOC in the sequel) on your toy problem.
  For a given value K of U(x, y) the level curves deduced from pdsolve or MOC solutions are the same, provided some care is taken (some level curves may not be real or a real level curve may have multiple branches).
  I'm not familiar with pdsolve but I used to use MOC several years ago. Depending on the initial conditions two characteristic curves (the level curves) may intersect, generating shocks where the solution has multiple values (one for each intersecting characteristic).
  I know how to handle these shocks with MOC but I don't know what kind of solution pdsolve builds: is it a regular solution or a generalized solution containing shocks?

Note that MOC does not contain any compatibility condition.
I have to think twice to the real reason which makes a compatibility condition to appear when I write

(4*x+2*y)*diff(u(x, y), x) + (x+3*y)*diff(u(x, y), y) = 0:

diff(u(x, y), x) = -(x+3*y);
diff(u(x, y), y) = (4*x+2*y);

as there are none when I write 

diff(x(s), s) = 4*x(s)+2*y(s);
diff(y(s), s) = x(s)+3*y(s);

I was taught of SA the first year after college, and of MOC a few years latter. I had forgotten the first one until I replied this question and mentioned SA because it was very simple. Until you pointed out this issue, I would have said that the two methods were equivalent for continuous solutions (ie shocks-free solutions).

@emendes 

Characteristic_curves.mw

See also Method_of_characteristics for a primer

Illustrative example

plots:-dualaxisplot(
  plot(x^2, x=-10..0, color=red, axis[2]=[color=white]), 
  plot(x^2, x=-10..0, color=red)
);

Explanation: with dualaxisplot the vertical axis of the second plot is always set at the right of this plot and its tickmarks are right oriented.
To make disappear the vertical axis of the first plot set its color to white.

workaround.mw

@Aixleft aimbot 

Ypu can get all the zeros of denom(v1) (and thus denom(v2)) this way

 

data := {alpha = 1/2 + 1/10*sqrt(5), beta = -1/2 + 1/10*sqrt(5)}:

# at this point get and tr still contain alpha and beta

det := simplify(N1*N4 - N2*N3, size):
tr := simplify(N1 + N4, size):

# To get the true value of det:

eval(det, data);

# at this point v1 still contains alpha and beta

v1 := simplify((N1+N4-sqrt((N1+N4)^2-4*(N1*N4-N2*N3)))*(1/2), size):

# synthetic expression of denom(v1):

factor(combine(denom(v1), trig));

                 24*beta*alpha*(sin(beta*v)+sin(v*(2*alpha+beta)))

# true expression of denom(v1)

eval(%, data):

# Zeros of denom(v1)

DenomZeros := solve(%=0, v, allsolutions);

     DenomZeros := -10*Pi*(2*_Z2+_B2) / ((2*_B2-3+sqrt(5))*(5+3*sqrt(5)))

# Properties of the protected names Maple uses

ind := convert((indets(DenomZeros) minus {Pi}), list);

            ind := [_B2, _Z2]
about~(ind);

Originally _B2, renamed _B2~:
  is assumed to be: OrProp(0,1)  # means _BE is either 0 or 1

Originally _Z2, renamed _Z2~:
  is assumed to be: integer

# some exact values of the zeros of denom(v1)

DZ1 := { seq(seq(eval(DenomZeros, {_B2=b2, _Z2=z2}), z2=-5..5), b2=0..1) }

# some numerical values

DZ1 := evalf( { seq(seq(eval(DenomZeros, {_B2=b2, _Z2=z2}), z2=-5..5), b2=0..1) } )=
{-35.12407365, -28.09925892, -23.87865850, -21.07444419, 

  -19.53708422, -15.19550995, -14.04962946, -10.85393568, 

  -7.024814730, -6.512361408, -2.170787136, 0., 2.170787136, 

  6.512361408, 7.024814730, 10.85393568, 14.04962946, 

  15.19550995, 19.53708422, 21.07444419, 28.09925892, 35.12407365

  }

# some numerical values within a given range (here -20..20)

DZ1 := convert( select(z -> (z > -20 and z < 20), DZ1), list);

 [-19.53708422, -15.19550995, -14.04962946, -10.85393568, 

   -7.024814730, -6.512361408, -2.170787136, 0., 2.170787136, 

   6.512361408, 7.024814730, 10.85393568, 14.04962946, 

   15.19550995, 19.53708422]

eigenvalues_2_mmcdara.mw  (plots not displayed to avoid uploading a big file)

@Aixleft aimbot 

v₁ is equal to f(v)/g(v).
Here is the graph of g(v) in the range -20..20
The points where g(v)=0 are
 
and as no exact value of them all can be obtained here is the lists of their numerical values in the range -20..20
(note that g(v)=-g(-v))

-19.53708421, -15.19550994, -14.04962946, -10.85393567, -7.024814731, -6.512361403, -2.170787134, 0, 2.170787134, 6.512361403, 7.024814731, 10.85393567, 14.04962946, 15.19550994, 19.53708421

Unless f(v) vanishes at these points and that the limits f(v)/g(v) exist, this means there are 17 subranges within -20..20 where f(v)/g(v) is continuous, and so, potentially 34 points where v₁= -1 or -1 = v₁= 1.

The second point is that f(v) is equal to

where the N_i are expressions in v. Thus, still potentially, the term under the square root can be negative in the range you consider which mean v₁ can be complex.
Is that what you wnt to account for when you write abs(v₁) < 1?

For instance, in the range 6.512361403..7.024814731 (the two bounds are consecutive zeros of g(v)), abs(v₁) , that is the square root of v₁*conjugate(v₁)) is positive and its minimum value, 1.4051226623472968, can by found this way

Optimization:-NLPSolve(abs(v1), v=6.512361403 .. 7.024814731);
 [1.4051226623472968, [v = 7.024812738168562]]

The first thing to do is to select ranges among 

SearchRanges := [seq(DZ1[i]..DZ1[i+1], i=1..NZ-1, 1)];
[
   -20 .. -19.53708421, 
   -19.53708421 .. -15.19550994, 
   -15.19550994 .. -14.04962946, 
   -14.04962946 .. -10.85393567, 
   -10.85393567 .. -7.024814731, 
  ​​​​​​​ -7.024814731 .. -6.512361403, 
   -6.512361403 .. -2.170787134, 
  ​​​​​​​ -2.170787134 .. 0, 
   0 .. 2.170787134, 
  ​​​​​​​ 2.170787134 .. 6.512361403, 
   6.512361403 .. 7.024814731, 
  ​​​​​​​ 7.024814731 .. 10.85393567, 
   10.85393567 .. 14.04962946, 
  ​​​​​​​ 14.04962946 .. 15.19550994, 
   15.19550994 .. 19.53708421, 
   19.53708421 .. 20
]

where the minimum value of abs(v₁)  is lower than 1:

for r in SearchRanges do   
opt := Optimization:-NLPSolve(abs(v1), v=r);   
printf("range = %20a, min(abs(v1)) = %1.8e opt[1]  at %a\n", evalf[6](r), opt[1], evalf[6](opt[2][1])) 
end do:

# In red are ranges where abs(v1) > 1

range =     -20. .. -19.5371, min(abs(v1)) = 3.11700611e-01 opt[1]  at v = -20.0000
range = -19.5371 .. -15.1955, min(abs(v1)) = 9.37714009e-10 opt[1]  at v = -17.1518
range = -15.1955 .. -14.0496, min(abs(v1)) = 5.44301800e-02 opt[1]  at v = -14.0497
range = -14.0496 .. -10.8539, min(abs(v1)) = 5.70391423e-01 opt[1]  at v = -13.4196
range = -10.8539 .. -7.02481, min(abs(v1)) = 1.04252918e-07 opt[1]  at v = -7.97600
range = -7.02481 .. -6.51236, min(abs(v1)) = 1.40555184e+00 opt[1]  at v = -7.02481
range = -6.51236 .. -2.17079, min(abs(v1)) = 4.47154513e-09 opt[1]  at v = -4.04371
range = -2.17079 .. 0.      , min(abs(v1)) = 1.35602167e-08 opt[1]  at v = -1.49814
range =       0. .. 2.17079 , min(abs(v1)) = 1.35602167e-08 opt[1]  at v = 1.49814
range =  2.17079 .. 6.51236 , min(abs(v1)) = 4.47154513e-09 opt[1]  at v = 4.04371
range =  6.51236 .. 7.02481 , min(abs(v1)) = 1.40512266e+00 opt[1]  at v = 7.02481
range =  7.02481 .. 10.8539 , min(abs(v1)) = 1.04252918e-07 opt[1]  at v = 7.97600
range =  10.8539 .. 14.0496 , min(abs(v1)) = 5.70391423e-01 opt[1]  at v = 13.4196
range =  14.0496 .. 15.1955 , min(abs(v1)) = 5.40266382e-02 opt[1]  at v = 14.0496
range =  15.1955 .. 19.5371 , min(abs(v1)) = 9.37714009e-10 opt[1]  at v = 17.1518
range =   19.5371 .. 20.    , min(abs(v1)) = 3.11700600e-01 opt[1]  at v = 20.

In any acceptable range but the two extremes there are potentially two points where abs(v₁) = 1. For instance, within the second range
 

​​​​​​​
 

The range within SearchRanges[2] where abs(v1) < 1 can be found this way

PT0 := Optimization:-NLPSolve(abs(v1), v=SearchRanges[2]);
PT1 := Optimization:-NLPSolve((abs(v1)-1)^2, v=SearchRanges[2]);
PT2 := Optimization:-NLPSolve((abs(v1)-1)^2, v=eval(v, PT1[2]) .. op(2, SearchRanges[2]));

    PT0 := [(9.377140094154938e-10, [v = -17.15175332553186]]
    PT1 := [(6.630751886116464e-20, [v = -18.62510516020815]]
    PT2 :=  [(0.144423193958391, [v = t(-16.430585488049168]]

# As abs(v1) is decreasing in -18.62510..-17.1517 and increasing in -17.1517..-16.43058 the
# range where abs(v1) < 1 is

           v = (-18.62510516020815, -16.430585488049168)

The approach is detailed in eigenvalues_mmcdara.mw and two more examples are given.
But do not lure yourself by thinking that Maple can do the job without help: for each of the 14 "blue ranges" listed above you will have to guide Maple to get the desired ranges.

@Carl Love @nm

Thanks to both of you

@C_R 

Great, I vote up

@dharr 

Impressive, I vote up.

@imparter 

Here are several ways do achieve what you want: proposals.mw

Download disp.mw

NOTE: I don't know why the tickmarks font has been changed by plotsetup  (just do display(< Q1 | Q2 >) after plotsetup(default) to see what the jpg file should contain) ?
I'm not even sure this issue can be fixed because Plotting Devices help page says "Notes: Font control is not available in this driver. The Maple JPEG driver is based in part on software provided by the Independent JPEG Group. JPEG is designed for photographs of natural scenes, which unlike plots, do not have sharp edges and abrupt color changes. GIF is a superior file format for plots; JPEG will make sharp edges blurry.'

Any attempt to modify the tickmarks (using for instance 

axis[1]=[tickmarks=[seq(i=sprintf("%s", convert(i, string)), i=-1..1)]]

) leads to the message Warning, ignoring axis information in _AXIS[n] structures

@Ali Hassani 

Thanks for the precisions.

So, if I correctly understand, may I say you have a differential equation which models the behaviour of some physical system; you know from some physical expertise that its solution shoud be of the form, let's say, 1-exp(-5*t); but the use of DTM to get a series representation of the exact solution if the polynomial f(t) in your intial question?

Plotting f(t) for t > 0 shows clearly that f does not behave as Physics suggest it would do.
Maybe you should provide us the differential equation to try and see how we can do to help you.

I look forward to reading from you

First point
Over which t-range do you want to resume your polynomial by a sum like alpha[0]+add(alpha[i]*exp(-beta[i]*t), i=1.. 5) ?

As you wrote "sum of exponential functions with a negative power" I assume you mean that the betas are strictly negative, am I right?
In which case the range over which to f is to be approximated should be a..b, with a >=0.
Am I still right?

Second point
Even with a range  a..b, with a >=0.there exists some additional constraints for all the betas to be strictly negative.
For instance f(t) must be convex everywhere in a..b.
Look here to the zeroes of diff(f(t), t$2):

fsolve(diff(f(t), t$2)=0, t=-100..100);
plot(abs(diff(f(t), t$2)), t=-0.8..1.8, axis[2]=[mode=log])
     -0.5543081877, 0.5294373439, 0.5525008020, 1.539900777

Then , for instance, all the betas are strictly negative if you consider the range -0.5543081877.. 0.5294373439, but some are positive while others are negative in the range -0.5543081877.. 0.6,
Given the rapid variation of f(t) for |t| > 1  the local discrepancies between f(t) and its "sum of exponential" approximations arround the zeroes of diff(f(t), t$2) are graphically invisible... but they do exist.

In the range t=0..1, f is correctly approximated by 

0.299348428151280e-1+1.18657264326334*10^(-7)*exp(17.3637775324664*t)

(adding more exponential terms improves marginally the quality of the fit ovr this range).
The above approximation (denoted fit1 below) can be obtained this way

restart

f := 0.020399949322360296902872908942 + 0.0261353198432118595103693714851*t^3 + 0.0240968505875842806805439681431*t^4 + 0.0148456155621193706595799212802*t^5 + 0.0239969764160351203722354728376*t^2 + 0.0204278458408370651586217048716*t - 0.00450853634927256388740864146173*t^6 - 0.0355389767483113696513996149731*t^7 - 0.0766669789661906882315038416910*t^8 - 0.120843030849135239578151569663*t^9 - 0.153280689906711146639066606024*t^10 - 0.150288711858517536713273977277*t^11 - 0.0808171080937786380164380347445*t^12 + 0.0872390654213369913348407061899*t^13 + 0.373992140377042586618283139889*t^14 + 0.766807288928470485618700282187*t^15 + 1.19339994493571167326973251788*t^16 + 1.49476369302534328383069681700*t^17 + 1.41015598591182237492637420929*t^18 + 0.593451797299651247527539427688*t^19 - 1.31434443870999971750661332301*t^20:
f := unapply(f, t):

# Assuming this t-range
domain := 0..1;

g := unapply(alpha[0]+sum(alpha[i]*exp(-beta[i]*t), i=1.. n) , (t, n)) :

N := 1;

X := [seq](domain, 0.01):
Y := f~(X):

fit1 := Statistics:-NonlinearFit(g(t, N), X, Y, t);
plot([f(t), fit1], t=domain, color=[blue, red], legend=[typeset('f'(t)), typeset('g'(t, N))]);

# RSS stands for Residual Sum of Squares and gives a measure of the 
# quality of the approximation of f(t) by g(t, 1).
# As the step [0.01] goes to 0, RSS converges to the L2 norm of f(t)-g(t, 1)
# on the interval 0..1

RSS_1 := Statistics:-NonlinearFit(g(t, N), X, Y, t, output=residualsumofsquares);

Look to file Methods.mw for more informations and to see how you could proceed for higher order approximations (you will see that Statistics:-NonlinearFit has some limitations and that Maple provides better alternatives).

As the above file is restricted to domain(s) of the form a..b with a >= 0, you will find in file Methods_different_domains.mw some hints about how to find approximations over ranges a..b <=0  and ranges a..b where a > 0 and b > 0.

But this require to be extremely careful when beta constraints are written.