Repository: sympy/sympy

Issue:
Use general rules in integral transforms code
```
for example, laplace_transform(f(t).diff(t), t, x) should give LaplaceTransform(f(t), t, x) - f(0)

This should actually be fairly easy: there is a "TODO" comment in integrals/transforms/transforms.py function IntegralTransform.doit where calls to methods that are overridden in the derived classes can be inserted.
```

Original issue for #5813: http://code.google.com/p/sympy/issues/detail?id=2714
Original author: https://code.google.com/u/104531927090589914088/
Original owner: https://code.google.com/u/104531927090589914088/


Discussion:
github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [907895ac]       [3ce40809]
-       4.68±0.2s         293±10ms     0.06  polygon.PolygonArbitraryPoint.time_bench01

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
hanspi42: 
> The easiest is probably to define a doit method in LaplaceTransform, run the code you want and then, if needed, call super().doit() (or whatever the correct syntax is).

> An option may be to check isinstance(self, LaplaceTransform), but I do not know if it works. (May be even easier if it does.)

Afer many experiments and thoughts, I have now switched to writing an independent `doit()` for the Laplace transform. The order "try directly, then split into sum" is not optimum for Laplace, I have changed it to "split into sum, then work on sum terms", and it already gives better results for many Laplace transforms.

My next steps are:
- time shift rule
- differentiation rule

For differentiation of an undefined function `f(t)`, I will need to return an unevaluated expression that means "the n-th derivative evaluated at t=0". I don't know how to best do this, does anyone have am idea?  Would returning `Derivative(f(t), (t, n)).subs(t,0)` be OK?
oscarbenjamin: > Would returning `Derivative(f(t), (t, n)).subs(t,0)` be OK?

Yes, that's basically the only way to represent a derivative evaluated at a point.
oscarbenjamin: I think this was previously marked as a draft. Is this ready for review now?
hanspi42: > I think this was previously marked as a draft. Is this ready for review now?

Thanks for asking, it is not ready yet. I will tell you here when I think it is at a stage where a review makes sense.
hanspi42: @oscarbenjamin: it is now ready for a first review. Do you think this way of writing a `doit()` is OK? If you say yes, I will next implement the derivatives of undefined functions.
oscarbenjamin: > Do you think this way of writing a `doit()` is OK?

I just looked quickly but I think it looks okay. I would prefer to implement the bulk (or all) of the logic in a separate function that is called by doit rather than in the doit function itself but otherwise this seems fine.
hanspi42: > 
> 
> > Do you think this way of writing a `doit()` is OK?
> 
> I just looked quickly but I think it looks okay. I would prefer to implement the bulk (or all) of the logic in a separate function that is called by doit rather than in the doit function itself but otherwise this seems fine.

OK, I have changed it. I will put all rules for `AppliedUndef` into the new `_apply_rules_undef` and rules for concrete functions into `_compute_transform`.
hanspi42: I have now also added a first concrete rule. From here on it should be smooth flying ;)
hanspi42: I am now dealing with `exp` expressions and have seen that many convergence regions were wrong up to now. I am correcting this, but do we have somebody who'll be able to review that part from a mathematical point of view?
hanspi42: I found a way to introduce checks etc. into the rules. Now it is much shorter.
hanspi42: Very annoying: when I do `LT(cos(t*(w - 1)), t, s)`, the expression that is handed down to `doit()` is unfortunately `cos(t*w - t)`, and this does not match. `a*t`.

Does anybody know how I could automatically apply `.collect(t)` to the arguments of all functions in an expression?
hanspi42: @oscargus: I have implemented your suggestion, thanks again.

@oscarbenjamin: I found a solution, I first tidy up the expression with a new `_deep_collect(f, t)` function by traversing the tree leaves first and applying `.collect()` each time the function in the tree node is `Add()`.

I implemented it such that every table entry now has a fifth element, a one-argument lambda function that is applied before pattern matching. Like this other prep functions could be used as well.
oscarbenjamin: This is looking nice and also much faster than master:
```python
In [3]: %time laplace_transform(exp(-3*t), t, z)
CPU times: user 12 ms, sys: 0 ns, total: 12 ms
Wall time: 12.5 ms
Out[3]: 
⎛  1            ⎞
⎜─────, -3, True⎟
⎝z + 3          ⎠
```
We should try to structure the code in such a way that there is a function that can be called that will attempt rule-based evaluation of the Laplace transform and does not fall back on the integration code. Then the `laplace_transform` function or `LaplaceTransform.doit` can call that function and handle falling back to integration. That way it becomes possible for the integration code itself to make use of the rule-based Laplace transform function when a suitable integral arrives.
hanspi42: 
> We should try to structure the code in such a way that there is a function that can be called that will attempt rule-based evaluation of the Laplace transform and does not fall back on the integration code. Then the `laplace_transform` function or `LaplaceTransform.doit` can call that function and handle falling back to integration. That way it becomes possible for the integration code itself to make use of the rule-based Laplace transform function when a suitable integral arrives.

OK, I will prepare for that by having all rules in the same place. Let me first complete the rules; then we can do this step. I guess I'll need some guidance for doing this ;)

hanspi42: Now `diff` is implemented, it looks beuatiful ;))))

![diff](https://user-images.githubusercontent.com/19270101/142177656-d087172a-326b-43a6-ba1a-f93fcc691901.png)


oscarbenjamin: Yes, that is nice. I wonder if the pretty printer can be improved to look more like the latex printer:
```python
In [2]: s = Symbol('s')

In [3]: laplace_transform(diff(f(t), (t, 4)), t, s)
Out[3]: 
                                                                  ⎛  2      ⎞│      ⎛  3      ⎞│   
 4                                 3         2 ⎛d       ⎞│        ⎜ d       ⎟│      ⎜ d       ⎟│   
s ⋅LaplaceTransform(f(t), t, s) - s ⋅f(0) - s ⋅⎜──(f(t))⎟│    - s⋅⎜───(f(t))⎟│    - ⎜───(f(t))⎟│   
                                               ⎝dt      ⎠│t=0     ⎜  2      ⎟│      ⎜  3      ⎟│   
                                                                  ⎝dt       ⎠│t=0   ⎝dt       ⎠│t=0
```

hanspi42: Would be nice to have a better pretty printer. Meanwhile I have written a small helper function to solve differential equations.

Check my test fiile.

[0003_all_asserts.zip](https://github.com/sympy/sympy/files/7554363/0003_all_asserts.zip)


The code lines

'''
s1 = LT(diff(f(t), (t, 3)), t, s)
print('\n-->',s1)
ic = [c0, c1, c2, c3, c4]
s2 = laplace_transform_replace_ic(s1, f, F, t, s, ic)
print('-->',s2)
print('-->',solve(s2, F(s)))
'''

give 

'''
--> s**3*LaplaceTransform(f(t), t, s) - s**2*f(0) - s*Subs(Derivative(f(t), t), t, 0) - Subs(Derivative(f(t), (t, 2)), t, 0)
--> -c0*s**2 - c1*s - c2 + s**3*F(s)
--> [(c0*s**2 + c1*s + c2)/s**3]
'''

Before I code all checks, tests, etc: is this a proper way to go?
oscargus: I think the "problem" with the pretty printer is related to `⋅` and `⎜` and the used font which aren't the same width as the rest of the characters. Which is a bad thing when the printer relies on a fixed-width font and counting characters...

This is with DejaVu Sans Mono (from my editor, hence the vertical line and so on):
![image](https://user-images.githubusercontent.com/8114497/142294298-3475c14c-baa5-4299-80b4-c10f74237d9b.png)

Btw, really nice PR @hanspi42 !

oscarbenjamin: > I think the "problem" with the pretty printer is related to `⋅` and `⎜` and the used font which aren't the same width as the rest of the characters.

You seem to be paying much more attention to some things than me but much less to others! To me the problem is that whereas in the latex version the Laplace transform of an undefined function comes out like `L_t[f(t)](s)`  the pretty print version is `LaplaceTransform(f(t), t, s)`. I think that some people get confused about whether SymPy has "computed" the Laplace transform because of this difference (https://github.com/sympy/sympy/issues/7219#issuecomment-802092986).

Also yes, much kudos to @hanspi42 for making headway on this longstanding issue!
hanspi42: This is going to give a lot to review, so I will guide you reviewers through the review process. I have re-checked all rules from the top to 3.14 and think they are implemented correctly now. The pdf below shows the first four fields of each rule converted automatically to LaTeX for better readability:

[0006_table.pdf](https://github.com/sympy/sympy/files/7574933/0006_table.pdf)

Could somebody check this against the respective tables on http://eqworld.ipmnet.ru/en/auxiliary/inttrans/ ?

oscarbenjamin: Is the special handling of `DiracDelta` in `_laplace_transform` still needed after the changes in this PR?
hanspi42: > Is the special handling of `DiracDelta` in `_laplace_transform` still needed after the changes in this PR?

No, it is not, but I will build it back only after I have implemented all rules. I tried to remove it once and too many things broke at once.

Other question: I would like to add time scaling to all rules instead of implementing it as an additional step. This has the disadvantage that the same kind of code is in multiple places, but the distinct advantage that it is visible in the rule list. One example is here:

https://github.com/sympy/sympy/blob/49fbd49b631a25a6ce785afdb775321c0a3decac/sympy/integrals/transforms.py#L1363-L1366

In many rules this is alrerady built in naturally, such as

https://github.com/sympy/sympy/blob/49fbd49b631a25a6ce785afdb775321c0a3decac/sympy/integrals/transforms.py#L1335-L1338

Is it OK for you if I just include it everyhere where it makes sense? The function `dco` given in each rule makes this very easy to match.
oscarbenjamin: > Is it OK for you if I just include it everyhere where it makes sense?

I think so although I'm not sure exactly what the two alternatives you are suggesting are. Do you just mean that you will generalise the rules from [1] rather than apply them literally? In other words this rule would be change to match `1/sqrt(a*t + b)`:
https://github.com/sympy/sympy/blob/49fbd49b631a25a6ce785afdb775321c0a3decac/sympy/integrals/transforms.py#L1222-L1224

[1] http://eqworld.ipmnet.ru/en/auxiliary/aux-inttrans.htm
hanspi42: > In other words this rule would be change to match `1/sqrt(a*t + b)`:
> 
> https://github.com/sympy/sympy/blob/49fbd49b631a25a6ce785afdb775321c0a3decac/sympy/integrals/transforms.py#L1222-L1224

Yes, it would then be

```
    (1/sqrt(b*t+a),
     sqrt(pi/s)*exp(a/b*s)*erfc(sqrt(a/b*s))/b,
     b>0, S.Zero, dco),
```
oscarbenjamin: That seems reasonable. I'm not sure if there are other situations where a general delay handler would still be needed perhaps for delays involving undefined functions or something.
hanspi42: > > Is the special handling of `DiracDelta` in `_laplace_transform` still needed after the changes in this PR?
> 
> No, it is not, but I will build it back only after I have implemented all rules. I tried to remove it once and too many things broke at once.

Correction: yes, it is. The present code includes the masking property and can transform the following:

```
LT(exp(-t)*(DiracDelta(t)+DiracDelta(t-42)), t, s)
```

I do not see how I could implement the masking property with rules.
hanspi42: I found something very annoying: Test case `test_issue_14692` uses `I` as a factor to `t`.  Now I am at the limit of my mathemathical knowledge. I don't know all the conditions for all the rules in [1] except that they surely work if coefficients are real. I think the only way to play it safe is to check whether they are real and only apply the rules if this is so.

The best way I found to check whether something is real that will only be evaluated after `.xreplace` is this:

https://github.com/sympy/sympy/blob/c8c1981f4849021706d9880ed0f1b43dc3cfd05f/sympy/integrals/transforms.py#L1183

applied like that:

https://github.com/sympy/sympy/blob/c8c1981f4849021706d9880ed0f1b43dc3cfd05f/sympy/integrals/transforms.py#L1218-L1221

It works with `test_issue_14692`, but I will also apply it to `b`, of course.

I think I will add this to all rules. Then expressions with non-real coefficients wil be solved as an integral, and it is later possible to relax conditions once we know they can be relaxed for some rules.

What do you think about it?
hanspi42: Aaargh, the above does not work.  The code on line 1436 of

https://github.com/sympy/sympy/blob/c8c1981f4849021706d9880ed0f1b43dc3cfd05f/sympy/integrals/transforms.py#L1433-L1438

triggers in this case when the replace operation is successful.

I have tried many things but am stuck: does anyone have an idea how I can provide a sympy expression in the list that is to be taken by my code and then 

1. replace wildcards by values
2. evaluate and check whether true

Whatever I tried up to now, my `isreal()` function is evaluated before the wildcards are replaced.

Does anyone have an idea how to do this?


hanspi42: I think I have a solution to make the new code work recursively such that general rules can also apply to concrete functions.

https://github.com/sympy/sympy/blob/58efc63ad93041ee593a06a142c5c6466b0cd3d9/sympy/integrals/transforms.py#L1547-L1559

@oscarbenjamin: this solves the "why do we need `if not any(func.has(t) for func in f.atoms(AppliedUndef))`" that you asked about two weeks ago.
oscargus: > Does anyone have an idea how to do this?

Not sure this is applicable here, but would it be possible to make a `is_extended_real` check and if `None` return a Piecewise like `Piecewise((result, Eq(re(a), a)), (LaplaceTransform(input), True))`? Or `Eq(im(a), 0)`.

I checked the Laplace transform for 1/(b*x+a) with Mathematica and that returns: 

(E^((a s)/b) (Gamma[0,(a s)/b]+Log[a/b]+Log[b/a]))/b

not sure if that is just a rewrite of Ei or not. For b = I it returns a value separated into a CI and SI, which, as far as I can tell, is just a rewrite. However, it is also commonly stated that Ei is just for real-values and E1 should be used for complex arguments. Or that may just be the integral definition that only works for real-valued arguments and Ei is defined for complex arguments in SymPy. 

Maybe @jksuom can shed some light on this?
oscarbenjamin: We can also just have the pattern rules written with Piecewise where needed. The Piecewise will evaluate if the value is real:
```python
In [1]: p = Piecewise((1, Contains(x, Reals)), (2, True))

In [2]: p
Out[2]: 
⎧1  for x ∈ ℝ
⎨            
⎩2  otherwise

In [3]: p.subs(x, 4)
Out[3]: 1

```
hanspi42: Very good ideas. Meanwhile, I could do miracles if I was able to use a new hint called `applyrules` in addition to 'needeval' etc.

However, I could not find out in where and in which file the selection of allowed hints can be found. I tried code like

```
                hints['applyrules'] = False
                return k*LaplaceTransform(ma1[g].func(t), t, s, **hints)
```

but LaplaceTransform bails out with 

```
  File "sympy\core\function.py", line 284, in __new__
    raise ValueError("Unknown options: %s" % options)

ValueError: Unknown options: {'applyrules': False}
```

so I guess I have to add this hint name to some list somewhere.
hanspi42: > I checked the Laplace transform for 1/(b*x+a) with Mathematica and that returns:
> 
> (E^((a s)/b) (Gamma[0,(a s)/b]+Log[a/b]+Log[b/a]))/b
> 
> not sure if that is just a rewrite of Ei or not.

Well, for non-zero arguments, Log[a/b]+Log[b/a] is zero, and Gamma(0,x) is Ei(x), so it is the same as we have ;)
hanspi42: > For b = I it returns a value separated into a CI and SI, which, as far as I can tell, is just a rewrite.

I am currently planning to let all rules just work on parameters that I am really sure are correct; all other cases will be treated as before.  It will always be possible to extend the rules later, but I first need to have something not too big but consistent, or it may blow my mind ;))
hanspi42: Just in case you are wondering what I did with the hints: I added a hint that prevents endless loops and makes it possible to apply rules without evaluating integrals

https://github.com/sympy/sympy/blob/904e1926f52232d881741636b2cba21210f62262/sympy/integrals/transforms.py#L1502-L1504

hanspi42: Dear reviewers,

I think I have finished the task and the code is ready for review.

After the review, while doing necessary corrections, I will also
- [ ] update the documentation a bit,
- [ ] add an example for solving differential equations to the documentation,
- [ ] write proper notes,
- [ ] squash the commits.

I have omitted all general rules containing integrals (1.15 ... 1.25) as I do not really see how this would help anybody much; if you think people need them and can actually use them efficiently, tell me. 

Parameter ranges for some rules may be more restrictive than necessary, but I would very much prefer to just take it as we have it (and later increase parameter ranges, when necessary, with comments from the field) such that I can turn to working on the rules for the inverse Laplace transform next.

hanspi42: By the way, if you execute the tests like so,

```
import sympy
from sympy.integrals.tests.test_transforms import test_laplace_transform
sympy.SYMPY_DEBUG=True

import sys
file_path = '0007_debug_tests_out.txt'
sys.stderr = open(file_path, "w")
test_laplace_transform()
```

then you will see in the output file which rules were triggered for which test cases.  Here the output for your convenience:

[0007_debug_tests_out.zip](https://github.com/sympy/sympy/files/7625616/0007_debug_tests_out.zip)


hanspi42: > I have checked it and have nothing to object against (on the opposite: a nice, huge, step!).

Thank you ;)

> Some minor things for the sake of it:
> 
> Are all the "removed" tests just moved elsewhere? I cannot seem to find the cosh*sinh-test, but I may just not see it.

There are no removed tests, they are just re-sorted.  `cosh*sinh`, e.g., is here:

https://github.com/sympy/sympy/blob/b5191b805eba30488a54cc0e46dbab66052a0429/sympy/integrals/tests/test_transforms.py#L549

> It would be nice to add a comment for the tests that corresponds to close issues as it right now is not that clear without checking the issues etc.

Good point.  Should these be comment lines next to the tests, or can it be a description under "Other comments" together with the release notes?

oscargus: > > Some minor things for the sake of it:
> > Are all the "removed" tests just moved elsewhere? I cannot seem to find the cosh*sinh-test, but I may just not see it.
> 
> There are no removed tests, they are just re-sorted. `cosh*sinh`, e.g., is here:

Great!

> > It would be nice to add a comment for the tests that corresponds to close issues as it right now is not that clear without checking the issues etc.
> 
> Good point. Should these be comment lines next to the tests, or can it be a description under "Other comments" together with the release notes?

I'd say lines in the tests (as often there are separate, explicit, tests for the issues). If nothing else for later developers to know that the test is "special" and shouldn't be removed in the first place.

hanspi42: 
> > > It would be nice to add a comment for the tests that corresponds to close issues as it right now is not that clear without checking the issues etc.
 
> > 
> > Good point. Should these be comment lines next to the tests, or can it be a description under "Other comments" together with the release notes?

> I'd say lines in the tests (as often there are separate, explicit, tests for the issues). If nothing else for later developers to know that the test is "special" and shouldn't be removed in the first place.

I will do this.
hanspi42: My to-do-list for the final submission:

- [ ] Re-write `_laplace_apply_rules` such that it uses a keyword-only argument instead of a hint;
- [ ] make sure that the `doit` argument is passed on properly in recursive calls;
- [ ] re-write `_laplace_apply_rules` such that it returns `None` instead of raising an error and adapt `_compute_transform` to cover this;
- [ ] convert the exp, sin, cos, sinh, cosh rules to high-level rules;
- [ ] un-duplicate code in `doit`;
- [ ] add comments to those test cases that test issues fixed with this pull request;
- [ ] update the documentation;
- [ ] add an example for solving differential equations to the documentation;
- [ ] write proper release notes;
- [ ] squash the commits.
hanspi42: For the reviewers' convenience: the built-in rules automatically converted to LaTeX:

[0006_table.pdf](https://github.com/sympy/sympy/files/7648681/0006_table.pdf)


hanspi42: Dear @oscarbenjamin and @oscargus, I think I have finished all tasks. Once you give me the green lights, I will squash the commits and then signal to @oscarbenjamin that a merge of this pull request is possible.
oscarbenjamin: I think this is almost there. The main thing is it would be better to split out the repetitive sections from `_laplace_apply_rules` into separate named functions so that they could have their own docstrings and be tested separately. Right now looking at that function it's hard to see what it's doing from a high-level.
hanspi42: > I think this is almost there. The main thing is it would be better to split out the repetitive sections from `_laplace_apply_rules` into separate named functions so that they could have their own docstrings and be tested separately. Right now looking at that function it's hard to see what it's doing from a high-level.

I totally agree and have found a way to do it. 
hanspi42: Final task, squash commits, still open ;)
asmeurer: I'm happy to see this finally being done. I have a few comments:

- We should obviously be doing this for every transform, not just Laplace. Can we set up the basic infrastructure for all the transforms, even if we don't have any rules for them yet? I think once we have the design down, adding rules will be the easy part, and it will also be easy for people to contribute additional rules.
- It isn't clear to me if the rules are applied to the exact expression provided by the user or if the input is pre-processed. I think it would be a good idea to at least split the input expression linearly, so `laplace_transform(a*f + b*g)` only pattern matches against `f` and `g`. Otherwise even a simple constant factor will break the pattern match (if this is actually already being done, then pay this comment no mind).
- Let's design this in a way that makes it possible for users to inject their own rules if they want to. Should be as simple as doing something like

```py
laplace_rules = []
def _laplace_build_rules():
    rules.append((...))
    rules.append((...))
```

- Just a note that ideally we would use matchpy for this, but that requires a lot of other refactorings first. 
- You could also consider modeling this after the `manualintegrate` module. The way that module is implemented allows things like extracting the steps used (this is how the integral steps part of sympygamma.com works). OTOH manualintegrate is a bit more complex than this because here we generally only expect one rule to match, and any matching rule to always "work", which isn't the case with manualintegrate.
hanspi42: > I'm happy to see this finally being done. I have a few comments:
> 
> We should obviously be doing this for every transform, not just Laplace. Can we set up the basic infrastructure for all the transforms, even if we don't have any rules for them yet? I think once we have the design down, adding rules will be the easy part, and it will also be easy for people to contribute additional rules.

I will be happy to do this, but I really would like you to let me to do it in this order:

1. Finish this pull request without the possibility to let users inject their own rules.
2. In a second pull request, re-write the inverse Laplace transform such that it uses rules to deal with rational fractions. Then it is possible to solve linear differential equations by Laplace transform, and it may be ready for deployment in sympy 1.10
3. Then I will have enough knowledge to be able to participate in a discussion about how a basic infrastructure for all integral transfroms should be designed. So we specify an architecture together.
4. We (or I) implement that architecture, and in that step, we make it possible for people to append their own rules, and we also tidy up the messy code in `IntegralTransform` that prevents you from seeing the answer to the next question ;)

>  It isn't clear to me if the rules are applied to the exact expression provided by the user or if the input is pre-processed. I think it would be a good idea to at least split the input expression linearly, so `laplace_transform(a*f + b*g)` only pattern matches against `f` and `g`. Otherwise even a simple constant factor will break the pattern match (if this is actually already being done, then pay this comment no mind).

It is like this already, but not easy to read. So I would like to defer the following points also to step 4 of my list above.

>  Let's design this in a way that makes it possible for users to inject their own rules if they want to. 

>  Just a note that ideally we would use matchpy for this, but that requires a lot of other refactorings first.

>  You could also consider modeling this after the `manualintegrate` module. ...

Is it OK for you if I proceed like this?

asmeurer: I'm fine with doing things in further PRs. It's not like anything here is public API, so it's OK if we end up completely refactoring it to something different later. 
hanspi42: Very good, thanks, so I will reply to the rest of @oscarbenjamin's comments with that plan in mind.
hanspi42: Once more I believe I have everything now ;))

If you give me green lights, I will squash the commits.
oscarbenjamin: > It's not like anything here is public API

The PR in its current state adds a new public function `laplace_transform_replace_ic`. I am not sure about the interface of this function and whether it is the best way to implement the feature expected by users so I would rather postpone that change to a future PR so it can be discussed in detail on its own.
hanspi42: > > It's not like anything here is public API
> 
> The PR in its current state adds a new public function `laplace_transform_replace_ic`. I am not sure about the interface of this function and whether it is the best way to implement the feature expected by users so I would rather postpone that change to a future PR so it can be discussed in detail on its own.

So would I delete it again, or make it an interrnal helper function for the time being and document it as such?
oscarbenjamin: > So would I delete it again, or make it an interrnal helper function for the time being and document it as such?

Is it actually used by any of the other code here? AFAICT it is only used in the tests so it's not really an internal helper. I think I would rather remove it from this PR entirely.
oscarbenjamin: > I think I would rather remove it from this PR entirely.

But we can open a new issue or a separate PR to add this right now and discuss the details there.

Patch:
diff --git a/sympy/integrals/transforms.py b/sympy/integrals/transforms.py
index 70aa43f509..a6fbaa8707 100644
--- a/sympy/integrals/transforms.py
+++ b/sympy/integrals/transforms.py
@@ -1,14 +1,14 @@
 """ Integral Transforms """
 from functools import reduce, wraps
 from itertools import repeat
-
-from sympy.core import S, pi
+from sympy.core import S, pi, I
 from sympy.core.add import Add
-from sympy.core.function import (AppliedUndef, count_ops, expand,
-                                 expand_complex, expand_mul, Function, Lambda)
+from sympy.core.function import (AppliedUndef, count_ops, Derivative, expand,
+                                 expand_complex, expand_mul, Function, Lambda,
+                                 WildFunction)
 from sympy.core.mul import Mul
 from sympy.core.numbers import igcd, ilcm
-from sympy.core.relational import _canonical, Ge, Gt, Lt, Unequality
+from sympy.core.relational import _canonical, Ge, Gt, Lt, Unequality, Eq
 from sympy.core.sorting import default_sort_key, ordered
 from sympy.core.symbol import Dummy, symbols, Wild
 from sympy.core.traversal import postorder_traversal
@@ -16,14 +16,15 @@
 from sympy.functions.elementary.complexes import (re, arg, Abs, polar_lift,
                                                   periodic_argument)
 from sympy.functions.elementary.exponential import exp, log, exp_polar
-from sympy.functions.elementary.hyperbolic import cosh, coth, sinh, tanh
+from sympy.functions.elementary.hyperbolic import cosh, coth, sinh, tanh, asinh
 from sympy.functions.elementary.integers import ceiling
 from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
 from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
-from sympy.functions.elementary.trigonometric import cos, cot, sin, tan
-from sympy.functions.special.bessel import besselj
+from sympy.functions.elementary.trigonometric import cos, cot, sin, tan, atan
+from sympy.functions.special.bessel import besseli, besselj, besselk, bessely
 from sympy.functions.special.delta_functions import DiracDelta, Heaviside
-from sympy.functions.special.gamma_functions import gamma
+from sympy.functions.special.error_functions import erf, erfc, Ei
+from sympy.functions.special.gamma_functions import digamma, gamma, lowergamma
 from sympy.functions.special.hyper import meijerg
 from sympy.integrals import integrate, Integral
 from sympy.integrals.meijerint import _dummy
@@ -39,6 +40,7 @@
 from sympy.solvers.inequalities import _solve_inequality
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.iterables import iterable
+from sympy.utilities.misc import debug
 
 
 ##########################################################################
@@ -125,6 +127,23 @@ def _collapse_extra(self, extra):
             raise IntegralTransformError(self.__class__.name, None, '')
         return cond
 
+    def _try_directly(self, **hints):
+        T = None
+        try_directly = not any(func.has(self.function_variable)
+                               for func in self.function.atoms(AppliedUndef))
+        if try_directly:
+            try:
+                T = self._compute_transform(self.function,
+                    self.function_variable, self.transform_variable, **hints)
+            except IntegralTransformError:
+                T = None
+
+        fn = self.function
+        if not fn.is_Add:
+            fn = expand_mul(fn)
+        return fn, T
+
+
     def doit(self, **hints):
         """
         Try to evaluate the transform in closed form.
@@ -147,18 +166,13 @@ def doit(self, **hints):
         ``(simplify, noconds, needeval) = (True, False, False)``.
         """
         needeval = hints.pop('needeval', False)
-        try_directly = not any(func.has(self.function_variable)
-                               for func in self.function.atoms(AppliedUndef))
-        if try_directly:
-            try:
-                return self._compute_transform(self.function,
-                    self.function_variable, self.transform_variable, **hints)
-            except IntegralTransformError:
-                pass
+        simplify = hints.pop('simplify', True)
+        hints['simplify'] = simplify
 
-        fn = self.function
-        if not fn.is_Add:
-            fn = expand_mul(fn)
+        fn, T = self._try_directly(**hints)
+
+        if T is not None:
+            return T
 
         if fn.is_Add:
             hints['needeval'] = needeval
@@ -176,7 +190,10 @@ def doit(self, **hints):
                 elif len(x) > 2:
                     # some region parameters and a condition (Mellin, Laplace)
                     extra += [x[1:]]
-            res = Add(*ress)
+            if simplify==True:
+                res = Add(*ress).simplify()
+            else:
+                res = Add(*ress)
             if not extra:
                 return res
             try:
@@ -840,6 +857,9 @@ def fundamental_strip(self):
         return a, b
 
     def _compute_transform(self, F, s, x, **hints):
+        # IntegralTransform's doit will cause this hint to exist, but
+        # InverseMellinTransform should ignore it
+        hints.pop('simplify', True)
         global _allowed
         if _allowed is None:
             _allowed = {
@@ -1014,9 +1034,14 @@ def expand_dirac_delta(expr):
     """
     return _lin_eq2dict(expr, expr.atoms(DiracDelta))
 
+
 @_noconds
 def _laplace_transform(f, t, s_, simplify=True):
-    """ The backend function for Laplace transforms. """
+    """ The backend function for Laplace transforms.
+
+    This backend assumes that the frontend has already split sums
+    such that `f` is to an addition anymore.
+    """
     s = Dummy('s')
     a = Wild('a', exclude=[t])
     deltazero = []
@@ -1137,6 +1162,624 @@ def sbs(expr):
         aux = _simplifyconds(aux, s, a)
     return _simplify(F.subs(s, s_), simplify), sbs(a), _canonical(sbs(aux))
 
+def _laplace_deep_collect(f, t):
+    """
+    This is an internal helper function that traverses through the epression
+    tree of `f(t)` and collects arguments. The purpose of it is that
+    anything like `f(w*t-1*t-c)` will be written as `f((w-1)*t-c)` such that
+    it can match `f(a*t+b)`.
+    """
+    func = f.func
+    args = list(f.args)
+    if len(f.args) == 0:
+        return f
+    else:
+        for k in range(len(args)):
+            args[k] = _laplace_deep_collect(args[k], t)
+        if func.is_Add:
+            return func(*args).collect(t)
+        else:
+            return func(*args)
+
+def _laplace_build_rules(t, s):
+    """
+    This is an internal helper function that returns the table of Laplace
+    transfrom rules in terms of the time variable `t` and the frequency
+    variable `s`.  It is used by `_laplace_apply_rules`.
+    """
+    a = Wild('a', exclude=[t])
+    b = Wild('b', exclude=[t])
+    n = Wild('n', exclude=[t])
+    tau = Wild('tau', exclude=[t])
+    omega = Wild('omega', exclude=[t])
+    dco = lambda f: _laplace_deep_collect(f,t)
+    laplace_transform_rules = [
+    # ( time domain,
+    #   laplace domain,
+    #   condition, convergence plane, preparation function )
+    #
+    # Catch constant (would otherwise be treated by 2.12)
+    (a, a/s, S.true, S.Zero, dco),
+    # DiracDelta rules
+    (DiracDelta(a*t-b),
+     exp(-s*b/a)/Abs(a),
+     Or(And(a>0, b>=0), And(a<0, b<=0)), S.Zero, dco),
+    (DiracDelta(a*t-b),
+     S(0),
+     Or(And(a<0, b>=0), And(a>0, b<=0)), S.Zero, dco),
+    # Rules from http://eqworld.ipmnet.ru/en/auxiliary/inttrans/
+    # 2.1
+    (1,
+     1/s,
+     S.true, S.Zero, dco),
+    # 2.2 expressed in terms of Heaviside
+    (Heaviside(a*t-b),
+     exp(-s*b/a)/s,
+     And(a>0, b>0), S.Zero, dco),
+    (Heaviside(a*t-b),
+     (1-exp(-s*b/a))/s,
+     And(a<0, b<0), S.Zero, dco),
+    (Heaviside(a*t-b),
+     1/s,
+     And(a>0, b<=0), S.Zero, dco),
+    (Heaviside(a*t-b),
+     0,
+     And(a<0, b>0), S.Zero, dco),
+    # 2.3
+    (t,
+     1/s**2,
+     S.true, S.Zero, dco),
+    # 2.4
+    (1/(a*t+b),
+     -exp(-b/a*s)*Ei(-b/a*s)/a,
+     a>0, S.Zero, dco),
+    # 2.5 and 2.6 are covered by 2.11
+    # 2.7
+    (1/sqrt(a*t+b),
+     sqrt(a*pi/s)*exp(b/a*s)*erfc(sqrt(b/a*s))/a,
+     a>0, S.Zero, dco),
+    # 2.8
+    (sqrt(t)/(t+b),
+     sqrt(pi/s)-pi*sqrt(b)*exp(b*s)*erfc(sqrt(b*s)),
+     S.true, S.Zero, dco),
+    # 2.9
+    ((a*t+b)**(-S(3)/2),
+     2*b**(-S(1)/2)-2*(pi*s/a)**(S(1)/2)*exp(b/a*s)*erfc(sqrt(b/a*s))/a,
+     a>0, S.Zero, dco),
+    # 2.10
+    (t**(S(1)/2)*(t+a)**(-1),
+     (pi/s)**(S(1)/2)-pi*a**(S(1)/2)*exp(a*s)*erfc(sqrt(a*s)),
+     S.true, S.Zero, dco),
+    # 2.11
+    (1/(a*sqrt(t) + t**(3/2)),
+     pi*a**(S(1)/2)*exp(a*s)*erfc(sqrt(a*s)),
+     S.true, S.Zero, dco),
+    # 2.12
+    (t**n,
+     gamma(n+1)/s**(n+1),
+     n>-1, S.Zero, dco),
+    # 2.13
+    ((a*t+b)**n,
+     lowergamma(n+1, b/a*s)*exp(-b/a*s)/s**(n+1)/a,
+     And(n>-1, a>0), S.Zero, dco),
+    # 2.14
+    (t**n/(t+a),
+     a**n*gamma(n+1)*lowergamma(-n,a*s),
+     n>-1, S.Zero, dco),
+    # 3.1
+    (exp(a*t-tau),
+     exp(-tau)/(s-a),
+     S.true, a, dco),
+    # 3.2
+    (t*exp(a*t-tau),
+     exp(-tau)/(s-a)**2,
+     S.true, a, dco),
+    # 3.3
+    (t**n*exp(a*t),
+     gamma(n+1)/(s-a)**(n+1),
+     n>-1, a, dco),
+    # 3.4 and 3.5 cannot be covered here because they are
+    #     sums and only the individual sum terms will get here.
+    # 3.6
+    (exp(-a*t**2),
+     sqrt(pi/4/a)*exp(s**2/4/a)*erfc(s/sqrt(4*a)),
+     a>0, S.Zero, dco),
+    # 3.7
+    (t*exp(-a*t**2),
+     1/(2*a)-2/sqrt(pi)/(4*a)**(S(3)/2)*s*erfc(s/sqrt(4*a)),
+     S.true, S.Zero, dco),
+    # 3.8
+    (exp(-a/t),
+     2*sqrt(a/s)*besselk(1, 2*sqrt(a*s)),
+     a>=0, S.Zero, dco),
+    # 3.9
+    (sqrt(t)*exp(-a/t),
+     S(1)/2*sqrt(pi/s**3)*(1+2*sqrt(a*s))*exp(-2*sqrt(a*s)),
+     a>=0, S.Zero, dco),
+    # 3.10
+    (exp(-a/t)/sqrt(t),
+     sqrt(pi/s)*exp(-2*sqrt(a*s)),
+     a>=0, S.Zero, dco),
+    # 3.11
+    (exp(-a/t)/(t*sqrt(t)),
+     sqrt(pi/a)*exp(-2*sqrt(a*s)),
+     a>0, S.Zero, dco),
+    # 3.12
+    (t**n*exp(-a/t),
+     2*(a/s)**((n+1)/2)*besselk(n+1, 2*sqrt(a*s)),
+     a>0, S.Zero, dco),
+    # 3.13
+    (exp(-2*sqrt(a*t)),
+     s**(-1)-sqrt(pi*a)*s**(-S(3)/2)*exp(a/s)*erfc(sqrt(a/s)),
+     S.true, S.Zero, dco),
+    # 3.14
+    (exp(-2*sqrt(a*t))/sqrt(t),
+     (pi/s)**(S(1)/2)*exp(a/s)*erfc(sqrt(a/s)),
+     S.true, S.Zero, dco),
+    # 4.1
+    (sinh(a*t),
+     a/(s**2-a**2),
+     S.true, Abs(a), dco),
+    # 4.2
+    (sinh(a*t)**2,
+     2*a**2/(s**3-4*a**2*s**2),
+     S.true, Abs(2*a), dco),
+    # 4.3
+    (sinh(a*t)/t,
+     log((s+a)/(s-a))/2,
+     S.true, a, dco),
+    # 4.4
+    (t**n*sinh(a*t),
+     gamma(n+1)/2*((s-a)**(-n-1)-(s+a)**(-n-1)),
+     n>-2, Abs(a), dco),
+    # 4.5
+    (sinh(2*sqrt(a*t)),
+     sqrt(pi*a)/s/sqrt(s)*exp(a/s),
+     S.true, S.Zero, dco),
+    # 4.6
+    (sqrt(t)*sinh(2*sqrt(a*t)),
+     pi**(S(1)/2)*s**(-S(5)/2)*(s/2+a)*exp(a/s)*erf(sqrt(a/s))-a**(S(1)/2)*s**(-2),
+     S.true, S.Zero, dco),
+    # 4.7
+    (sinh(2*sqrt(a*t))/sqrt(t),
+     pi**(S(1)/2)*s**(-S(1)/2)*exp(a/s)*erf(sqrt(a/s)),
+     S.true, S.Zero, dco),
+    # 4.8
+    (sinh(sqrt(a*t))**2/sqrt(t),
+     pi**(S(1)/2)/2*s**(-S(1)/2)*(exp(a/s)-1),
+     S.true, S.Zero, dco),
+    # 4.9
+    (cosh(a*t),
+     s/(s**2-a**2),
+     S.true, Abs(a), dco),
+    # 4.10
+    (cosh(a*t)**2,
+     (s**2-2*a**2)/(s**3-4*a**2*s**2),
+     S.true, Abs(2*a), dco),
+    # 4.11
+    (t**n*cosh(a*t),
+     gamma(n+1)/2*((s-a)**(-n-1)+(s+a)**(-n-1)),
+     n>-1, Abs(a), dco),
+    # 4.12
+    (cosh(2*sqrt(a*t)),
+     1/s+sqrt(pi*a)/s/sqrt(s)*exp(a/s)*erf(sqrt(a/s)),
+     S.true, S.Zero, dco),
+    # 4.13
+    (sqrt(t)*cosh(2*sqrt(a*t)),
+     pi**(S(1)/2)*s**(-S(5)/2)*(s/2+a)*exp(a/s),
+     S.true, S.Zero, dco),
+    # 4.14
+    (cosh(2*sqrt(a*t))/sqrt(t),
+     pi**(S(1)/2)*s**(-S(1)/2)*exp(a/s),
+     S.true, S.Zero, dco),
+    # 4.15
+    (cosh(sqrt(a*t))**2/sqrt(t),
+     pi**(S(1)/2)/2*s**(-S(1)/2)*(exp(a/s)+1),
+     S.true, S.Zero, dco),
+    # 5.1
+    (log(a*t),
+     -log(s/a+S.EulerGamma)/s,
+     a>0, S.Zero, dco),
+    # 5.2
+    (log(1+a*t),
+     -exp(s/a)/s*Ei(-s/a),
+     S.true, S.Zero, dco),
+    # 5.3
+    (log(a*t+b),
+     (log(b)-exp(s/b/a)/s*a*Ei(-s/b))/s*a,
+     a>0, S.Zero, dco),
+    # 5.4 is covered by 5.7
+    # 5.5
+    (log(t)/sqrt(t),
+     -sqrt(pi/s)*(log(4*s)+S.EulerGamma),
+     S.true, S.Zero, dco),
+    # 5.6 is covered by 5.7
+    # 5.7
+    (t**n*log(t),
+     gamma(n+1)*s**(-n-1)*(digamma(n+1)-log(s)),
+     n>-1, S.Zero, dco),
+    # 5.8
+    (log(a*t)**2,
+     ((log(s/a)+S.EulerGamma)**2+pi**2/6)/s,
+     a>0, S.Zero, dco),
+    # 5.9
+    (exp(-a*t)*log(t),
+     -(log(s+a)+S.EulerGamma)/(s+a),
+     S.true, -a, dco),
+    # 6.1
+    (sin(omega*t),
+     omega/(s**2+omega**2),
+     S.true, S.Zero, dco),
+    # 6.2
+    (Abs(sin(omega*t)),
+     omega/(s**2+omega**2)*coth(pi*s/2/omega),
+     omega>0, S.Zero, dco),
+    # 6.3 and 6.4 are covered by 1.8
+    # 6.5 is covered by 1.8 together with 2.5
+    # 6.6
+    (sin(omega*t)/t,
+     atan(omega/s),
+     S.true, S.Zero, dco),
+    # 6.7
+    (sin(omega*t)**2/t,
+     log(1+4*omega**2/s**2)/4,
+     S.true, S.Zero, dco),
+    # 6.8
+    (sin(omega*t)**2/t**2,
+     omega*atan(2*omega/s)-s*log(1+4*omega**2/s**2)/4,
+     S.true, S.Zero, dco),
+    # 6.9
+    (sin(2*sqrt(a*t)),
+      sqrt(pi*a)/s/sqrt(s)*exp(-a/s),
+      a>0, S.Zero, dco),
+    # 6.10
+    (sin(2*sqrt(a*t))/t,
+     pi*erf(sqrt(a/s)),
+     a>0, S.Zero, dco),
+    # 6.11
+    (cos(omega*t),
+     s/(s**2+omega**2),
+     S.true, S.Zero, dco),
+    # 6.12
+    (cos(omega*t)**2,
+     (s**2+2*omega**2)/(s**2+4*omega**2)/s,
+     S.true, S.Zero, dco),
+    # 6.13 is covered by 1.9 together with 2.5
+    # 6.14 and 6.15 cannot be done with this method, the respective sum
+    #       parts do not converge. Solve elsewhere if really needed.
+    # 6.16
+    (sqrt(t)*cos(2*sqrt(a*t)),
+     sqrt(pi)/2*s**(-S(5)/2)*(s-2*a)*exp(-a/s),
+     a>0, S.Zero, dco),
+    # 6.17
+    (cos(2*sqrt(a*t))/sqrt(t),
+     sqrt(pi/s)*exp(-a/s),
+     a>0, S.Zero, dco),
+    # 6.18
+    (sin(a*t)*sin(b*t),
+     2*a*b*s/(s**2+(a+b)**2)/(s**2+(a-b)**2),
+     S.true, S.Zero, dco),
+    # 6.19
+    (cos(a*t)*sin(b*t),
+     b*(s**2-a**2+b**2)/(s**2+(a+b)**2)/(s**2+(a-b)**2),
+     S.true, S.Zero, dco),
+    # 6.20
+    (cos(a*t)*cos(b*t),
+     s*(s**2+a**2+b**2)/(s**2+(a+b)**2)/(s**2+(a-b)**2),
+     S.true, S.Zero, dco),
+    # 6.21
+    (exp(b*t)*sin(a*t),
+     a/((s-b)**2+a**2),
+     S.true, b, dco),
+    # 6.22
+    (exp(b*t)*cos(a*t),
+     (s-b)/((s-b)**2+a**2),
+     S.true, b, dco),
+    # 7.1
+    (erf(a*t),
+     exp(s**2/(2*a)**2)*erfc(s/(2*a))/s,
+     a>0, S.Zero, dco),
+    # 7.2
+    (erf(sqrt(a*t)),
+     sqrt(a)/sqrt(s+a)/s,
+     a>0, S.Zero, dco),
+    # 7.3
+    (exp(a*t)*erf(sqrt(a*t)),
+     sqrt(a)/sqrt(s)/(s-a),
+     a>0, a, dco),
+    # 7.4
+    (erf(sqrt(a/t)/2),
+     (1-exp(-sqrt(a*s)))/s,
+     a>0, S.Zero, dco),
+    # 7.5
+    (erfc(sqrt(a*t)),
+     (sqrt(s+a)-sqrt(a))/sqrt(s+a)/s,
+     a>0, S.Zero, dco),
+    # 7.6
+    (exp(a*t)*erfc(sqrt(a*t)),
+     1/(s+sqrt(a*s)),
+     a>0, S.Zero, dco),
+    # 7.7
+    (erfc(sqrt(a/t)/2),
+     exp(-sqrt(a*s))/s,
+     a>0, S.Zero, dco),
+    # 8.1, 8.2
+    (besselj(n, a*t),
+     a**n/(sqrt(s**2+a**2)*(s+sqrt(s**2+a**2))**n),
+     And(a>0, n>-1), S.Zero, dco),
+    # 8.3, 8.4
+    (t**b*besselj(n, a*t),
+     2**n/sqrt(pi)*gamma(n+S.Half)*a**n*(s**2+a**2)**(-n-S.Half),
+     And(And(a>0, n>-S.Half), Eq(b, n)), S.Zero, dco),
+    # 8.5
+    (t**b*besselj(n, a*t),
+     2**(n+1)/sqrt(pi)*gamma(n+S(3)/2)*a**n*s*(s**2+a**2)**(-n-S(3)/2),
+     And(And(a>0, n>-1), Eq(b, n+1)), S.Zero, dco),
+    # 8.6
+    (besselj(0, 2*sqrt(a*t)),
+     exp(-a/s)/s,
+     a>0, S.Zero, dco),
+    # 8.7, 8.8
+    (t**(b)*besselj(n, 2*sqrt(a*t)),
+     a**(n/2)*s**(-n-1)*exp(-a/s),
+     And(And(a>0, n>-1), Eq(b, n*S.Half)), S.Zero, dco),
+    # 8.9
+    (besselj(0, a*sqrt(t**2+b*t)),
+     exp(b*s-b*sqrt(s**2+a**2))/sqrt(s**2+a**2),
+     b>0, S.Zero, dco),
+    # 8.10, 8.11
+    (besseli(n, a*t),
+     a**n/(sqrt(s**2-a**2)*(s+sqrt(s**2-a**2))**n),
+     And(a>0, n>-1), Abs(a), dco),
+    # 8.12
+    (t**b*besseli(n, a*t),
+     2**n/sqrt(pi)*gamma(n+S.Half)*a**n*(s**2-a**2)**(-n-S.Half),
+     And(And(a>0, n>-S.Half), Eq(b, n)), Abs(a), dco),
+    # 8.13
+    (t**b*besseli(n, a*t),
+     2**(n+1)/sqrt(pi)*gamma(n+S(3)/2)*a**n*s*(s**2-a**2)**(-n-S(3)/2),
+     And(And(a>0, n>-1), Eq(b, n+1)), Abs(a), dco),
+    # 8.15, 8.16
+    (t**(b)*besseli(n, 2*sqrt(a*t)),
+     a**(n/2)*s**(-n-1)*exp(a/s),
+     And(And(a>0, n>-1), Eq(b, n*S.Half)), S.Zero, dco),
+    # 8.17
+    (bessely(0, a*t),
+     -2/pi*asinh(s/a)/sqrt(s**2+a**2),
+     a>0, S.Zero, dco),
+    # 8.18
+    (besselk(0, a*t),
+     (log(s+sqrt(s**2-a**2)))/(sqrt(s**2-a**2)),
+     a>0, Abs(a), dco)
+    ]
+    return laplace_transform_rules
+
+def _laplace_cr(f, a, c, **hints):
+    """
+    Internal helper function that will return `(f, a, c)` unless `**hints`
+    contains `noconds=True`, in which case it will only return `f`.
+    """
+    conds = not hints.get('noconds', False)
+    if conds:
+        return f, a, c
+    else:
+        return f
+
+def _laplace_rule_timescale(f, t, s, doit=True, **hints):
+    r"""
+    This internal helper function tries to apply the time-scaling rule of the
+    Laplace transform and returns `None` if it cannot do it.
+
+    Time-scaling means the following: if $F(s)$ is the Laplace transform of,
+    $f(t)$, then, for any $a>0$, the Laplace transform of $f(at)$ will be
+    $\frac1a F(\frac{s}{a})$. This scaling will also affect the transform's
+    convergence plane.
+    """
+    _simplify = hints.pop('simplify', True)
+    b = Wild('b', exclude=[t])
+    g = WildFunction('g', nargs=1)
+    k, func = f.as_independent(t, as_Add=False)
+    ma1 = func.match(g)
+    if ma1:
+        arg = ma1[g].args[0].collect(t)
+        ma2 = arg.match(b*t)
+        if ma2 and ma2[b]>0:
+            debug('_laplace_apply_rules match:')
+            debug('      f:    %s ( %s, %s )'%(f, ma1, ma2))
+            debug('      rule: amplitude and time scaling (1.1, 1.2)')
+            if ma2[b]==1:
+                if doit==True and not any(func.has(t) for func
+                                          in ma1[g].atoms(AppliedUndef)):
+                    return k*_laplace_transform(ma1[g].func(t), t, s,
+                                                simplify=_simplify)
+                else:
+                    return k*LaplaceTransform(ma1[g].func(t), t, s, **hints)
+            else:
+                L = _laplace_apply_rules(ma1[g].func(t), t, s/ma2[b],
+                                         doit=doit, **hints)
+                try:
+                    r, p, c = L
+                    return (k/ma2[b]*r, p, c)
+                except TypeError:
+                    return k/ma2[b]*L
+    return None
+
+def _laplace_rule_heaviside(f, t, s, doit=True, **hints):
+    """
+    This internal helper function tries to transform a product containing the
+    `Heaviside` function and returns `None` if it cannot do it.
+    """
+    hints.pop('simplify', True)
+    a = Wild('a', exclude=[t])
+    b = Wild('b', exclude=[t])
+    y = Wild('y')
+    g = WildFunction('g', nargs=1)
+    k, func = f.as_independent(t, as_Add=False)
+    ma1 = func.match(Heaviside(y)*g)
+    if ma1:
+        ma2 = ma1[y].match(t-a)
+        ma3 = ma1[g].args[0].collect(t).match(t-b)
+        if ma2 and ma2[a]>0 and ma3 and ma2[a]==ma3[b]:
+            debug('_laplace_apply_rules match:')
+            debug('      f:    %s ( %s, %s, %s )'%(f, ma1, ma2, ma3))
+            debug('      rule: time shift (1.3)')
+            L = _laplace_apply_rules(ma1[g].func(t), t, s, doit=doit, **hints)
+            try:
+                r, p, c = L
+                return (k*exp(-ma2[a]*s)*r, p, c)
+            except TypeError:
+                return k*exp(-ma2[a]*s)*L
+    return None
+
+
+def _laplace_rule_exp(f, t, s, doit=True, **hints):
+    """
+    This internal helper function tries to transform a product containing the
+    `exp` function and returns `None` if it cannot do it.
+    """
+    hints.pop('simplify', True)
+    a = Wild('a', exclude=[t])
+
+    y = Wild('y')
+    z = Wild('z')
+    k, func = f.as_independent(t, as_Add=False)
+    ma1 = func.match(exp(y)*z)
+    if ma1:
+        ma2 = ma1[y].collect(t).match(a*t)
+        if ma2:
+            debug('_laplace_apply_rules match:')
+            debug('      f:    %s ( %s, %s )'%(f, ma1, ma2))
+            debug('      rule: multiply with exp (1.5)')
+            L = _laplace_apply_rules(ma1[z], t, s-ma2[a], doit=doit, **hints)
+            try:
+                r, p, c = L
+                return (r, p+ma2[a], c)
+            except TypeError:
+                return L
+    return None
+
+def _laplace_rule_trig(f, t, s, doit=True, **hints):
+    """
+    This internal helper function tries to transform a product containing a
+    trigonometric function (`sin`, `cos`, `sinh`, `cosh`, ) and returns
+    `None` if it cannot do it.
+    """
+    _simplify = hints.pop('simplify', True)
+    a = Wild('a', exclude=[t])
+    y = Wild('y')
+    z = Wild('z')
+    k, func = f.as_independent(t, as_Add=False)
+    # All of the rules have a very similar form: trig(y)*z is matched, and then
+    # two copies of the Laplace transform of z are shifted in the s Domain
+    # and added with a weight; see rules 1.6 to 1.9 in
+    # http://eqworld.ipmnet.ru/en/auxiliary/inttrans/laplace1.pdf
+    # The parameters in the tuples are (fm, nu, s1, s2, sd):
+    #   fm: Function to match
+    #   nu: Number of the rule, for debug purposes
+    #   s1: weight of the sum, 'I' for sin and '1' for all others
+    #   s2: sign of the second copy of the Laplace transform of z
+    #   sd: shift direction; shift along real or imaginary axis if `1` or `I`
+    trigrules = [(sinh(y), '1.6',  1, -1, 1), (cosh(y), '1.7', 1, 1, 1),
+                 (sin(y),  '1.8', -I, -1, I), (cos(y),  '1.9', 1, 1, I)]
+    for trigrule in trigrules:
+        fm, nu, s1, s2, sd = trigrule
+        ma1 = func.match(fm*z)
+        if ma1:
+            ma2 = ma1[y].collect(t).match(a*t)
+            if ma2:
+                debug('_laplace_apply_rules match:')
+                debug('      f:    %s ( %s, %s )'%(f, ma1, ma2))
+                debug('      rule: multiply with %s (%s)'%(fm.func, nu))
+                L = _laplace_apply_rules(ma1[z], t, s, doit=doit, **hints)
+                try:
+                    r, p, c = L
+                    # The convergence plane changes only if the shift has been
+                    # done along the real axis:
+                    if sd==1:
+                        cp_shift = Abs(ma2[a])
+                    else:
+                        cp_shift = 0
+                    return ((s1*(r.subs(s, s-sd*ma2[a])+\
+                                    s2*r.subs(s, s+sd*ma2[a]))).simplify()/2,
+                            p+cp_shift, c)
+                except TypeError:
+                    if doit==True and _simplify==True:
+                        return (s1*(L.subs(s, s-sd*ma2[a])+\
+                                    s2*L.subs(s, s+sd*ma2[a]))).simplify()/2
+                    else:
+                        return (s1*(L.subs(s, s-sd*ma2[a])+\
+                                    s2*L.subs(s, s+sd*ma2[a])))/2
+    return None
+
+def _laplace_rule_diff(f, t, s, doit=True, **hints):
+    """
+    This internal helper function tries to transform an expression containing
+    a derivative of an undefined function and returns `None` if it cannot
+    do it.
+    """
+    hints.pop('simplify', True)
+    a = Wild('a', exclude=[t])
+    y = Wild('y')
+    n = Wild('n', exclude=[t])
+    g = WildFunction('g', nargs=1)
+    ma1 = f.match(a*Derivative(g, (t, n)))
+    if ma1 and ma1[g].args[0] == t and ma1[n].is_integer:
+        debug('_laplace_apply_rules match:')
+        debug('      f:    %s'%(f,))
+        debug('      rule: time derivative (1.11, 1.12)')
+        d = []
+        for k in range(ma1[n]):
+            if k==0:
+                y = ma1[g].func(t).subs(t, 0)
+            else:
+                y = Derivative(ma1[g].func(t), (t, k)).subs(t, 0)
+            d.append(s**(ma1[n]-k-1)*y)
+        r = s**ma1[n]*_laplace_apply_rules(ma1[g].func(t), t, s, doit=doit,
+                                           **hints)
+        return r - Add(*d)
+    return None
+
+
+def _laplace_apply_rules(f, t, s, doit=True, **hints):
+    """
+    Helper function for the class LaplaceTransform.
+
+    This function does a Laplace transform based on rules and, after
+    applying the rules, hands the rest over to `_laplace_transform`, which
+    will attempt to integrate.
+
+    If it is called with `doit=False`, then it will instead return
+    `LaplaceTransform` objects.
+    """
+
+    k, func = f.as_independent(t, as_Add=False)
+
+    simple_rules = _laplace_build_rules(t, s)
+    for t_dom, s_dom, check, plane, prep in simple_rules:
+        ma = prep(func).match(t_dom)
+        if ma:
+            debug('_laplace_apply_rules match:')
+            debug('      f:    %s'%(func,))
+            debug('      rule: %s o---o %s'%(t_dom, s_dom))
+            try:
+                debug('      try   %s'%(check,))
+                c = check.xreplace(ma)
+                debug('      check %s -> %s'%(check, c))
+                if c==True:
+                    return _laplace_cr(k*s_dom.xreplace(ma),
+                                plane.xreplace(ma), S.true, **hints)
+            except Exception:
+                debug('_laplace_apply_rules did not match.')
+    if f.has(DiracDelta):
+        return None
+
+    prog_rules = [_laplace_rule_timescale, _laplace_rule_heaviside,
+                  _laplace_rule_exp, _laplace_rule_trig, _laplace_rule_diff]
+    for p_rule in prog_rules:
+        LT = p_rule(f, t, s, doit=doit, **hints)
+        if LT is not None:
+            return LT
+    return None
 
 class LaplaceTransform(IntegralTransform):
     """
@@ -1151,7 +1794,14 @@ class LaplaceTransform(IntegralTransform):
     _name = 'Laplace'
 
     def _compute_transform(self, f, t, s, **hints):
-        return _laplace_transform(f, t, s, **hints)
+        LT = _laplace_apply_rules(f, t, s, **hints)
+        if LT is None:
+            _simplify = hints.pop('simplify', True)
+            debug('_laplace_apply_rules could not match function %s'%(f,))
+            debug('    hints: %s'%(hints,))
+            return _laplace_transform(f, t, s, simplify=_simplify, **hints)
+        else:
+            return LT
 
     def _as_integral(self, f, t, s):
         return Integral(f*exp(-s*t), (t, S.Zero, S.Infinity))
@@ -1169,6 +1819,20 @@ def _collapse_extra(self, extra):
                 'Laplace', None, 'No combined convergence.')
         return plane, cond
 
+    def _try_directly(self, **hints):
+        fn = self.function
+        debug('----> _try_directly: %s'%(fn, ))
+        t_ = self.function_variable
+        s_ = self.transform_variable
+        LT = None
+        if not fn.is_Add:
+            fn = expand_mul(fn)
+            try:
+                LT = self._compute_transform(fn, t_, s_, **hints)
+            except IntegralTransformError:
+                LT = None
+        return fn, LT
+
 
 def laplace_transform(f, t, s, legacy_matrix=True, **hints):
     r"""
@@ -1180,21 +1844,29 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints):
     ===========
 
     For all sensible functions, this converges absolutely in a
-    half plane  `a < \operatorname{Re}(s)`.
+    half-plane
+
+    .. math :: a < \operatorname{Re}(s)
 
     This function returns ``(F, a, cond)`` where ``F`` is the Laplace
-    transform of ``f``, `\operatorname{Re}(s) > a` is the half-plane
-    of convergence, and ``cond`` are auxiliary convergence conditions.
+    transform of ``f``, `a` is the half-plane of convergence, and `cond` are
+    auxiliary convergence conditions.
+
+    The implementation is rule-based, and if you are interested in which
+    rules are applied, and whether integration is attemped, you can switch
+    debug information on by setting `sympy.SYMPY_DEBUG=True`.
 
-    The lower bound is `0^{-}`, meaning that this bound should be approached
+    The lower bound is `0-`, meaning that this bound should be approached
     from the lower side. This is only necessary if distributions are involved.
     At present, it is only done if `f(t)` contains ``DiracDelta``, in which
-    case the Laplace transform is computed as
+    case the Laplace transform is computed implicitly as
 
-    .. math :: F(s) = \lim_{\tau\to 0^{-}} \int_{\tau}^\infty e^{-st} f(t) \mathrm{d}t.
+    .. math :: F(s) = \lim_{\tau\to 0^{-}} \int_{\tau}^\infty e^{-st} f(t) \mathrm{d}t
 
-    If the integral cannot be computed in closed form, this function returns
-    an unevaluated :class:`LaplaceTransform` object.
+    by applying rules.
+
+    If the integral cannot be fully computed in closed form, this function
+    returns an unevaluated :class:`LaplaceTransform` object.
 
     For a description of possible hints, refer to the docstring of
     :func:`sympy.integrals.transforms.IntegralTransform.doit`. If ``noconds=True``,
@@ -1211,13 +1883,13 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints):
     Examples
     ========
 
-    >>> from sympy.integrals import laplace_transform
+    >>> from sympy.integrals.transforms import laplace_transform
     >>> from sympy.abc import t, s, a
     >>> from sympy.functions import DiracDelta, exp
-    >>> laplace_transform(t**a, t, s)
-    (gamma(a + 1)/(s*s**a), 0, re(a) > -1)
+    >>> laplace_transform(t**4, t, s)
+    (24/s**5, 0, True)
     >>> laplace_transform(DiracDelta(t)-a*exp(-a*t),t,s)
-    (-a/(a + s) + 1, 0, Abs(arg(a)) <= pi/2)
+    (s/(a + s), Max(0, -a), True)
 
     See Also
     ========
@@ -1227,6 +1899,8 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints):
 
     """
 
+    debug('\n***** laplace_transform(%s, %s, %s)'%(f, t, s))
+
     if isinstance(f, MatrixBase) and hasattr(f, 'applyfunc'):
 
         conds = not hints.get('noconds', False)


Test patch:
diff --git a/sympy/integrals/tests/test_transforms.py b/sympy/integrals/tests/test_transforms.py
index 79413f1316..1b89e5ddac 100644
--- a/sympy/integrals/tests/test_transforms.py
+++ b/sympy/integrals/tests/test_transforms.py
@@ -1,6 +1,6 @@
 from sympy.integrals.transforms import (mellin_transform,
-    inverse_mellin_transform, laplace_transform, inverse_laplace_transform,
-    fourier_transform, inverse_fourier_transform,
+    inverse_mellin_transform, laplace_transform,
+    inverse_laplace_transform, fourier_transform, inverse_fourier_transform,
     sine_transform, inverse_sine_transform,
     cosine_transform, inverse_cosine_transform,
     hankel_transform, inverse_hankel_transform,
@@ -8,20 +8,20 @@
     InverseLaplaceTransform, InverseFourierTransform,
     InverseSineTransform, InverseCosineTransform, IntegralTransformError)
 from sympy.core.function import (Function, expand_mul)
-from sympy.core import EulerGamma
+from sympy.core import EulerGamma, Subs, Derivative, diff
 from sympy.core.numbers import (I, Rational, oo, pi)
-from sympy.core.relational import Eq, Ne
+from sympy.core.relational import Eq
 from sympy.core.singleton import S
 from sympy.core.symbol import (Symbol, symbols)
 from sympy.functions.combinatorial.factorials import factorial
-from sympy.functions.elementary.complexes import (Abs, arg, re, unpolarify)
+from sympy.functions.elementary.complexes import (Abs, re, unpolarify)
 from sympy.functions.elementary.exponential import (exp, exp_polar, log)
-from sympy.functions.elementary.hyperbolic import (cosh, sinh)
+from sympy.functions.elementary.hyperbolic import (cosh, sinh, coth, asinh)
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (atan, atan2, cos, sin, tan)
 from sympy.functions.special.bessel import (besseli, besselj, besselk, bessely)
 from sympy.functions.special.delta_functions import Heaviside
-from sympy.functions.special.error_functions import (erf, erfc, expint)
+from sympy.functions.special.error_functions import (erf, erfc, expint, Ei)
 from sympy.functions.special.gamma_functions import gamma
 from sympy.functions.special.hyper import meijerg
 from sympy.simplify.gammasimp import gammasimp
@@ -38,10 +38,7 @@ def test_undefined_function():
     f = Function('f')
     assert mellin_transform(f(x), x, s) == MellinTransform(f(x), x, s)
     assert mellin_transform(f(x) + exp(-x), x, s) == \
-        (MellinTransform(f(x), x, s) + gamma(s), (0, oo), True)
-
-    assert laplace_transform(2*f(x), x, s) == 2*LaplaceTransform(f(x), x, s)
-    # TODO test derivative and other rules when implemented
+        (MellinTransform(f(x), x, s) + gamma(s + 1)/s, (0, oo), True)
 
 
 def test_free_symbols():
@@ -468,71 +465,263 @@ def mysimp(expr):
 
 @slow
 def test_laplace_transform():
+    from sympy import lowergamma
     from sympy.functions.special.delta_functions import DiracDelta
     from sympy.functions.special.error_functions import (fresnelc, fresnels)
     LT = laplace_transform
-    a, b, c, = symbols('a b c', positive=True)
-    t = symbols('t')
-    w = Symbol("w")
+    a, b, c, = symbols('a, b, c', positive=True)
+    t, w, x = symbols('t, w, x')
     f = Function("f")
-
-    # Test unevaluated form
-    assert laplace_transform(f(t), t, w) == LaplaceTransform(f(t), t, w)
+    g = Function("g")
+
+    # Test rule-base evaluation according to
+    # http://eqworld.ipmnet.ru/en/auxiliary/inttrans/
+    # Power-law functions (laplace2.pdf)
+    assert LT(a*t+t**2+t**(S(5)/2), t, s) ==\
+        (a/s**2 + 2/s**3 + 15*sqrt(pi)/(8*s**(S(7)/2)), 0, True)
+    assert LT(b/(t+a), t, s) == (-b*exp(-a*s)*Ei(-a*s), 0, True)
+    assert LT(1/sqrt(t+a), t, s) ==\
+        (sqrt(pi)*sqrt(1/s)*exp(a*s)*erfc(sqrt(a)*sqrt(s)), 0, True)
+    assert LT(sqrt(t)/(t+a), t, s) ==\
+        (-pi*sqrt(a)*exp(a*s)*erfc(sqrt(a)*sqrt(s)) + sqrt(pi)*sqrt(1/s),
+         0, True)
+    assert LT((t+a)**(-S(3)/2), t, s) ==\
+        (-2*sqrt(pi)*sqrt(s)*exp(a*s)*erfc(sqrt(a)*sqrt(s)) + 2/sqrt(a),
+         0, True)
+    assert LT(t**(S(1)/2)*(t+a)**(-1), t, s) ==\
+        (-pi*sqrt(a)*exp(a*s)*erfc(sqrt(a)*sqrt(s)) + sqrt(pi)*sqrt(1/s),
+         0, True)
+    assert LT(1/(a*sqrt(t) + t**(3/2)), t, s) ==\
+        (pi*sqrt(a)*exp(a*s)*erfc(sqrt(a)*sqrt(s)), 0, True)
+    assert LT((t+a)**b, t, s) ==\
+        (s**(-b - 1)*exp(-a*s)*lowergamma(b + 1, a*s), 0, True)
+    assert LT(t**5/(t+a), t, s) == (120*a**5*lowergamma(-5, a*s), 0, True)
+    # Exponential functions (laplace3.pdf)
+    assert LT(exp(t), t, s) == (1/(s - 1), 1, True)
+    assert LT(exp(2*t), t, s) == (1/(s - 2), 2, True)
+    assert LT(exp(a*t), t, s) == (1/(s - a), a, True)
+    assert LT(exp(a*(t-b)), t, s) == (exp(-a*b)/(-a + s), a, True)
+    assert LT(t*exp(-a*(t)), t, s) == ((a + s)**(-2), -a, True)
+    assert LT(t*exp(-a*(t-b)), t, s) == (exp(a*b)/(a + s)**2, -a, True)
+    assert LT(b*t*exp(-a*t), t, s) == (b/(a + s)**2, -a, True)
+    assert LT(t**(S(7)/4)*exp(-8*t)/gamma(S(11)/4), t, s) ==\
+        ((s + 8)**(-S(11)/4), -8, True)
+    assert LT(t**(S(3)/2)*exp(-8*t), t, s) ==\
+        (3*sqrt(pi)/(4*(s + 8)**(S(5)/2)), -8, True)
+    assert LT(t**a*exp(-a*t), t, s) ==  ((a+s)**(-a-1)*gamma(a+1), -a, True)
+    assert LT(b*exp(-a*t**2), t, s) ==\
+        (sqrt(pi)*b*exp(s**2/(4*a))*erfc(s/(2*sqrt(a)))/(2*sqrt(a)), 0, True)
+    assert LT(exp(-2*t**2), t, s) ==\
+        (sqrt(2)*sqrt(pi)*exp(s**2/8)*erfc(sqrt(2)*s/4)/4, 0, True)
+    assert LT(b*exp(2*t**2), t, s) == b*LaplaceTransform(exp(2*t**2), t, s)
+    assert LT(t*exp(-a*t**2), t, s) ==\
+        (1/(2*a) - s*erfc(s/(2*sqrt(a)))/(4*sqrt(pi)*a**(S(3)/2)), 0, True)
+    assert LT(exp(-a/t), t, s) ==\
+        (2*sqrt(a)*sqrt(1/s)*besselk(1, 2*sqrt(a)*sqrt(s)), 0, True)
+    assert LT(sqrt(t)*exp(-a/t), t, s) ==\
+        (sqrt(pi)*(2*sqrt(a)*sqrt(s) + 1)*sqrt(s**(-3))*exp(-2*sqrt(a)*\
+                                                    sqrt(s))/2, 0, True)
+    assert LT(exp(-a/t)/sqrt(t), t, s) ==\
+        (sqrt(pi)*sqrt(1/s)*exp(-2*sqrt(a)*sqrt(s)), 0, True)
+    assert LT( exp(-a/t)/(t*sqrt(t)), t, s) ==\
+        (sqrt(pi)*sqrt(1/a)*exp(-2*sqrt(a)*sqrt(s)), 0, True)
+    assert LT(exp(-2*sqrt(a*t)), t, s) ==\
+        ( 1/s -sqrt(pi)*sqrt(a) * exp(a/s)*erfc(sqrt(a)*sqrt(1/s))/\
+         s**(S(3)/2), 0, True)
+    assert LT(exp(-2*sqrt(a*t))/sqrt(t), t, s) == (exp(a/s)*erfc(sqrt(a)*\
+        sqrt(1/s))*(sqrt(pi)*sqrt(1/s)), 0, True)
+    assert LT(t**4*exp(-2/t), t, s) ==\
+        (8*sqrt(2)*(1/s)**(S(5)/2)*besselk(5, 2*sqrt(2)*sqrt(s)), 0, True)
+    # Hyperbolic functions (laplace4.pdf)
+    assert LT(sinh(a*t), t, s) == (a/(-a**2 + s**2), a, True)
+    assert LT(b*sinh(a*t)**2, t, s) == (2*a**2*b/(-4*a**2*s**2 + s**3),
+                                        2*a, True)
+    # The following line confirms that issue #21202 is solved
+    assert LT(cosh(2*t), t, s) == (s/(-4 + s**2), 2, True)
+    assert LT(cosh(a*t), t, s) == (s/(-a**2 + s**2), a, True)
+    assert LT(cosh(a*t)**2, t, s) == ((-2*a**2 + s**2)/(-4*a**2*s**2 + s**3),
+                                      2*a, True)
+    assert LT(sinh(x + 3), x, s) == (
+        (-s + (s + 1)*exp(6) + 1)*exp(-3)/(s - 1)/(s + 1)/2, 0, Abs(s) > 1)
+    # The following line replaces the old test test_issue_7173()
+    assert LT(sinh(a*t)*cosh(a*t), t, s) == (a/(-4*a**2 + s**2), 2*a, True)
+    assert LT(sinh(a*t)/t, t, s) == (log((a + s)/(-a + s))/2, a, True)
+    assert LT(t**(-S(3)/2)*sinh(a*t), t, s) ==\
+        (-sqrt(pi)*(sqrt(-a + s) - sqrt(a + s)), a, True)
+    assert LT(sinh(2*sqrt(a*t)), t, s) ==\
+        (sqrt(pi)*sqrt(a)*exp(a/s)/s**(S(3)/2), 0, True)
+    assert LT(sqrt(t)*sinh(2*sqrt(a*t)), t, s) ==\
+        (-sqrt(a)/s**2 + sqrt(pi)*(a + s/2)*exp(a/s)*erf(sqrt(a)*\
+                                            sqrt(1/s))/s**(S(5)/2), 0, True)
+    assert LT(sinh(2*sqrt(a*t))/sqrt(t), t, s) ==\
+        (sqrt(pi)*exp(a/s)*erf(sqrt(a)*sqrt(1/s))/sqrt(s), 0, True)
+    assert LT(sinh(sqrt(a*t))**2/sqrt(t), t, s) ==\
+        (sqrt(pi)*(exp(a/s) - 1)/(2*sqrt(s)), 0, True)
+    assert LT(t**(S(3)/7)*cosh(a*t), t, s) ==\
+        (((a + s)**(-S(10)/7) + (-a+s)**(-S(10)/7))*gamma(S(10)/7)/2, a, True)
+    assert LT(cosh(2*sqrt(a*t)), t, s) ==\
+        (sqrt(pi)*sqrt(a)*exp(a/s)*erf(sqrt(a)*sqrt(1/s))/s**(S(3)/2) + 1/s,
+         0, True)
+    assert LT(sqrt(t)*cosh(2*sqrt(a*t)), t, s) ==\
+        (sqrt(pi)*(a + s/2)*exp(a/s)/s**(S(5)/2), 0, True)
+    assert LT(cosh(2*sqrt(a*t))/sqrt(t), t, s) ==\
+        (sqrt(pi)*exp(a/s)/sqrt(s), 0, True)
+    assert LT(cosh(sqrt(a*t))**2/sqrt(t), t, s) ==\
+        (sqrt(pi)*(exp(a/s) + 1)/(2*sqrt(s)), 0, True)
+    # logarithmic functions (laplace5.pdf)
+    assert LT(log(t), t, s) == (-log(s+S.EulerGamma)/s, 0, True)
+    assert LT(log(t/a), t, s) == (-log(a*s + S.EulerGamma)/s, 0, True)
+    assert LT(log(1+a*t), t, s) == (-exp(s/a)*Ei(-s/a)/s, 0, True)
+    assert LT(log(t+a), t, s) == ((log(a) - exp(s/a)*Ei(-s/a)/s)/s, 0, True)
+    assert LT(log(t)/sqrt(t), t, s) ==\
+        (sqrt(pi)*(-log(s) - 2*log(2) - S.EulerGamma)/sqrt(s), 0, True)
+    assert LT(t**(S(5)/2)*log(t), t, s) ==\
+        (15*sqrt(pi)*(-log(s)-2*log(2)-S.EulerGamma+S(46)/15)/(8*s**(S(7)/2)),
+         0, True)
+    assert (LT(t**3*log(t), t, s, noconds=True)-6*(-log(s) - S.EulerGamma\
+                                    + S(11)/6)/s**4).simplify() == S.Zero
+    assert LT(log(t)**2, t, s) ==\
+        (((log(s) + EulerGamma)**2 + pi**2/6)/s, 0, True)
+    assert LT(exp(-a*t)*log(t), t, s) ==\
+        ((-log(a + s) - S.EulerGamma)/(a + s), -a, True)
+    # Trigonometric functions (laplace6.pdf)
+    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
+    assert LT(Abs(sin(a*t)), t, s) ==\
+        (a*coth(pi*s/(2*a))/(a**2 + s**2), 0, True)
+    assert LT(sin(a*t)/t, t, s) == (atan(a/s), 0, True)
+    assert LT(sin(a*t)**2/t, t, s) == (log(4*a**2/s**2 + 1)/4, 0, True)
+    assert LT(sin(a*t)**2/t**2, t, s) ==\
+        (a*atan(2*a/s) - s*log(4*a**2/s**2 + 1)/4, 0, True)
+    assert LT(sin(2*sqrt(a*t)), t, s) ==\
+        (sqrt(pi)*sqrt(a)*exp(-a/s)/s**(S(3)/2), 0, True)
+    assert LT(sin(2*sqrt(a*t))/t, t, s) == (pi*erf(sqrt(a)*sqrt(1/s)), 0, True)
+    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
+    assert LT(cos(a*t)**2, t, s) ==\
+        ((2*a**2 + s**2)/(s*(4*a**2 + s**2)), 0, True)
+    assert LT(sqrt(t)*cos(2*sqrt(a*t)), t, s) ==\
+        (sqrt(pi)*(-2*a + s)*exp(-a/s)/(2*s**(S(5)/2)), 0, True)
+    assert LT(cos(2*sqrt(a*t))/sqrt(t), t, s) ==\
+        (sqrt(pi)*sqrt(1/s)*exp(-a/s), 0, True)
+    assert LT(sin(a*t)*sin(b*t), t, s) ==\
+        (2*a*b*s/((s**2 + (a - b)**2)*(s**2 + (a + b)**2)), 0, True)
+    assert LT(cos(a*t)*sin(b*t), t, s) ==\
+        (b*(-a**2 + b**2 + s**2)/((s**2 + (a - b)**2)*(s**2 + (a + b)**2)),
+         0, True)
+    assert LT(cos(a*t)*cos(b*t), t, s) ==\
+        (s*(a**2 + b**2 + s**2)/((s**2 + (a - b)**2)*(s**2 + (a + b)**2)),
+         0, True)
+    assert LT(c*exp(-b*t)*sin(a*t), t, s) == (a*c/(a**2 + (b + s)**2),
+                                              -b, True)
+    assert LT(c*exp(-b*t)*cos(a*t), t, s) == ((b + s)*c/(a**2 + (b + s)**2),
+                                              -b, True)
+    assert LT(cos(x + 3), x, s) == ((s*cos(3) - sin(3))/(s**2 + 1), 0, True)
+    # Error functions (laplace7.pdf)
+    assert LT(erf(a*t), t, s) == (exp(s**2/(4*a**2))*erfc(s/(2*a))/s, 0, True)
+    assert LT(erf(sqrt(a*t)), t, s) == (sqrt(a)/(s*sqrt(a + s)), 0, True)
+    assert LT(exp(a*t)*erf(sqrt(a*t)), t, s) ==\
+        (sqrt(a)/(sqrt(s)*(-a + s)), a, True)
+    assert LT(erf(sqrt(a/t)/2), t, s) == ((1-exp(-sqrt(a)*sqrt(s)))/s, 0, True)
+    assert LT(erfc(sqrt(a*t)), t, s) ==\
+        ((-sqrt(a) + sqrt(a + s))/(s*sqrt(a + s)), 0, True)
+    assert LT(exp(a*t)*erfc(sqrt(a*t)), t, s) ==\
+        (1/(sqrt(a)*sqrt(s) + s), 0, True)
+    assert LT(erfc(sqrt(a/t)/2), t, s) == (exp(-sqrt(a)*sqrt(s))/s, 0, True)
+    # Bessel functions (laplace8.pdf)
+    assert LT(besselj(0, a*t), t, s) == (1/sqrt(a**2 + s**2), 0, True)
+    assert LT(besselj(1, a*t), t, s) ==\
+        (a/(sqrt(a**2 + s**2)*(s + sqrt(a**2 + s**2))), 0, True)
+    assert LT(besselj(2, a*t), t, s) ==\
+        (a**2/(sqrt(a**2 + s**2)*(s + sqrt(a**2 + s**2))**2), 0, True)
+    assert LT(t*besselj(0, a*t), t, s) ==\
+        (s/(a**2 + s**2)**(S(3)/2), 0, True)
+    assert LT(t*besselj(1, a*t), t, s) ==\
+        (a/(a**2 + s**2)**(S(3)/2), 0, True)
+    assert LT(t**2*besselj(2, a*t), t, s) ==\
+        (3*a**2/(a**2 + s**2)**(S(5)/2), 0, True)
+    assert LT(besselj(0, 2*sqrt(a*t)), t, s) == (exp(-a/s)/s, 0, True)
+    assert LT(t**(S(3)/2)*besselj(3, 2*sqrt(a*t)), t, s) ==\
+        (a**(S(3)/2)*exp(-a/s)/s**4, 0, True)
+    assert LT(besselj(0, a*sqrt(t**2+b*t)), t, s) ==\
+        (exp(b*s - b*sqrt(a**2 + s**2))/sqrt(a**2 + s**2), 0, True)
+    assert LT(besseli(0, a*t), t, s) == (1/sqrt(-a**2 + s**2), a, True)
+    assert LT(besseli(1, a*t), t, s) ==\
+        (a/(sqrt(-a**2 + s**2)*(s + sqrt(-a**2 + s**2))), a, True)
+    assert LT(besseli(2, a*t), t, s) ==\
+        (a**2/(sqrt(-a**2 + s**2)*(s + sqrt(-a**2 + s**2))**2), a, True)
+    assert LT(t*besseli(0, a*t), t, s) == (s/(-a**2 + s**2)**(S(3)/2), a, True)
+    assert LT(t*besseli(1, a*t), t, s) == (a/(-a**2 + s**2)**(S(3)/2), a, True)
+    assert LT(t**2*besseli(2, a*t), t, s) ==\
+        (3*a**2/(-a**2 + s**2)**(S(5)/2), a, True)
+    assert LT(t**(S(3)/2)*besseli(3, 2*sqrt(a*t)), t, s) ==\
+        (a**(S(3)/2)*exp(a/s)/s**4, 0, True)
+    assert LT(bessely(0, a*t), t, s) ==\
+        (-2*asinh(s/a)/(pi*sqrt(a**2 + s**2)), 0, True)
+    assert LT(besselk(0, a*t), t, s) ==\
+        (log(s + sqrt(-a**2 + s**2))/sqrt(-a**2 + s**2), a, True)
+    assert LT(sin(a*t)**8, t, s) ==\
+        (40320*a**8/(s*(147456*a**8 + 52480*a**6*s**2 + 4368*a**4*s**4 +\
+                        120*a**2*s**6 + s**8)), 0, True)
+
+    # Test general rules and unevaluated forms
+    # These all also test whether issue #7219 is solved.
+    assert LT(Heaviside(t-1)*cos(t-1), t, s) == (s*exp(-s)/(s**2 + 1), 0, True)
+    assert LT(a*f(t), t, w) == a*LaplaceTransform(f(t), t, w)
+    assert LT(a*Heaviside(t+1)*f(t+1), t, s) ==\
+        a*LaplaceTransform(f(t + 1)*Heaviside(t + 1), t, s)
+    assert LT(a*Heaviside(t-1)*f(t-1), t, s) ==\
+        a*LaplaceTransform(f(t), t, s)*exp(-s)
+    assert LT(b*f(t/a), t, s) == a*b*LaplaceTransform(f(t), t, a*s)
+    assert LT(exp(-f(x)*t), t, s) == (1/(s + f(x)), -f(x), True)
+    assert LT(exp(-a*t)*f(t), t, s) == LaplaceTransform(f(t), t, a + s)
+    assert LT(exp(-a*t)*erfc(sqrt(b/t)/2), t, s) ==\
+        (exp(-sqrt(b)*sqrt(a + s))/(a + s), -a, True)
+    assert LT(sinh(a*t)*f(t), t, s) ==\
+        LaplaceTransform(f(t), t, -a+s)/2 - LaplaceTransform(f(t), t, a+s)/2
+    assert LT(sinh(a*t)*t, t, s) ==\
+        (-1/(2*(a + s)**2) + 1/(2*(-a + s)**2), a, True)
+    assert LT(cosh(a*t)*f(t), t, s) ==\
+        LaplaceTransform(f(t), t, -a+s)/2 + LaplaceTransform(f(t), t, a+s)/2
+    assert LT(cosh(a*t)*t, t, s) ==\
+        (1/(2*(a + s)**2) + 1/(2*(-a + s)**2), a, True)
+    assert LT(sin(a*t)*f(t), t, s) ==\
+        I*(-LaplaceTransform(f(t), t, -I*a + s) +\
+           LaplaceTransform(f(t), t, I*a + s))/2
+    assert LT(sin(a*t)*t, t, s) ==\
+        (2*a*s/(a**4 + 2*a**2*s**2 + s**4), 0, True)
+    assert LT(cos(a*t)*f(t), t, s) ==\
+        LaplaceTransform(f(t), t, -I*a + s)/2 +\
+        LaplaceTransform(f(t), t, I*a + s)/2
+    assert LT(cos(a*t)*t, t, s) ==\
+        ((-a**2 + s**2)/(a**4 + 2*a**2*s**2 + s**4), 0, True)
+    # The following two lines test whether issues #5813 and #7176 are solved.
+    assert LT(diff(f(t), (t, 1)), t, s) == s*LaplaceTransform(f(t), t, s)\
+        - f(0)
+    assert LT(diff(f(t), (t, 3)), t, s) == s**3*LaplaceTransform(f(t), t, s)\
+        - s**2*f(0) - s*Subs(Derivative(f(t), t), t, 0)\
+            - Subs(Derivative(f(t), (t, 2)), t, 0)
+    assert LT(a*f(b*t)+g(c*t), t, s) == a*LaplaceTransform(f(t), t, s/b)/b +\
+        LaplaceTransform(g(t), t, s/c)/c
     assert inverse_laplace_transform(
         f(w), w, t, plane=0) == InverseLaplaceTransform(f(w), w, t, 0)
+    assert LT(f(t)*g(t), t, s) == LaplaceTransform(f(t)*g(t), t, s)
 
-    # test a bug
-    spos = symbols('s', positive=True)
-    assert LT(exp(t), t, spos) == (1/(spos - 1), 0, spos > 1)
-
-    # basic tests from wikipedia
+    # additional basic tests from wikipedia
     assert LT((t - a)**b*exp(-c*(t - a))*Heaviside(t - a), t, s) == \
         ((s + c)**(-b - 1)*exp(-a*s)*gamma(b + 1), -c, True)
-    assert LT(t**a, t, s) == (s**(-a - 1)*gamma(a + 1), 0, True)
-    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
-    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
-    assert LT(1 - exp(-a*t), t, s) == (a/(s*(a + s)), 0, True)
-
     assert LT((exp(2*t) - 1)*exp(-b - t)*Heaviside(t)/2, t, s, noconds=True) \
         == exp(-b)/(s**2 - 1)
 
-    assert LT(exp(t), t, s) == (1/(s - 1), 0, abs(s) > 1)
-    assert LT(exp(2*t), t, s) == (1/(s - 2), 0, abs(s) > 2)
-    assert LT(exp(a*t), t, s) == (1/(s - a), a, Ne(s/a, 1))
-
-    assert LT(log(t/a), t, s) == (-(log(a*s) + EulerGamma)/s, 0, True)
-
-    assert LT(erf(t), t, s) == (erfc(s/2)*exp(s**2/4)/s, 0, True)
-
-    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
-    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
-    # TODO would be nice to have these come out better
-    assert LT(exp(-a*t)*sin(b*t), t, s) == (b/(b**2 + (a + s)**2), -a, True)
-    assert LT(exp(-a*t)*cos(b*t), t, s) == \
-        ((a + s)/(b**2 + (a + s)**2), -a, True)
-
-    assert LT(besselj(0, t), t, s) == (1/sqrt(1 + s**2), 0, True)
-    assert LT(besselj(1, t), t, s) == (1 - 1/sqrt(1 + 1/s**2), 0, True)
-    # TODO general order works, but is a *mess*
-    # TODO besseli also works, but is an even greater mess
-
-    # test a bug in conditions processing
-    # TODO the auxiliary condition should be recognised/simplified
-    assert LT(exp(t)*cos(t), t, s)[:-1] in [
-        ((s - 1)/(s**2 - 2*s + 2), -oo),
-        ((s - 1)/((s - 1)**2 + 1), -oo),
-    ]
-
     # DiracDelta function: standard cases
-    assert LT(DiracDelta(t), t, s) == (1, -oo, True)
-    assert LT(DiracDelta(a*t), t, s) == (1/a, -oo, True)
-    assert LT(DiracDelta(t/42), t, s) == (42, -oo, True)
-    assert LT(DiracDelta(t+42), t, s) == (0, -oo, True)
+    assert LT(DiracDelta(t), t, s) == (1, 0, True)
+    assert LT(DiracDelta(a*t), t, s) == (1/a, 0, True)
+    assert LT(DiracDelta(t/42), t, s) == (42, 0, True)
+    assert LT(DiracDelta(t+42), t, s) == (0, 0, True)
     assert LT(DiracDelta(t)+DiracDelta(t-42), t, s) == \
-        (1 + exp(-42*s), -oo, True)
-    assert LT(DiracDelta(t)-a*exp(-a*t), t, s) == (-a/(a + s) + 1, 0, True)
+        (1 + exp(-42*s), 0, True)
+    assert LT(DiracDelta(t)-a*exp(-a*t), t, s) == (s/(a + s), 0, True)
     assert LT(exp(-t)*(DiracDelta(t)+DiracDelta(t-42)), t, s) == \
         (exp(-42*s - 42) + 1, -oo, True)
+
     # Collection of cases that cannot be fully evaluated and/or would catch
     # some common implementation errors
     assert LT(DiracDelta(t**2), t, s) == LaplaceTransform(DiracDelta(t**2), t, s)
@@ -542,8 +731,17 @@ def test_laplace_transform():
     assert LT((DiracDelta(t) + 1)*(DiracDelta(t - 1) + 1), t, s) == \
         (LaplaceTransform(DiracDelta(t)*DiracDelta(t - 1), t, s) + \
          1 + exp(-s) + 1/s, 0, True)
-    assert LT(DiracDelta(2*t - 2*exp(a)), t, s) == \
-        (exp(-s*exp(a))/2, -oo, True)
+    assert LT(DiracDelta(2*t-2*exp(a)), t, s) == (exp(-s*exp(a))/2, 0, True)
+    assert LT(DiracDelta(-2*t+2*exp(a)), t, s) == (exp(-s*exp(a))/2, 0, True)
+
+    # Heaviside tests
+    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
+    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
+    assert LT(Heaviside(t-1), t, s) == (exp(-s)/s, 0, True)
+    assert LT(Heaviside(2*t-4), t, s) == (exp(-2*s)/s, 0, True)
+    assert LT(Heaviside(-2*t+4), t, s) == ((1 - exp(-2*s))/s, 0, True)
+    assert LT(Heaviside(2*t+4), t, s) == (1/s, 0, True)
+    assert LT(Heaviside(-2*t+4), t, s) == ((1 - exp(-2*s))/s, 0, True)
 
     # Fresnel functions
     assert laplace_transform(fresnels(t), t, s) == \
@@ -553,6 +751,7 @@ def test_laplace_transform():
         ((2*sin(s**2/(2*pi))*fresnelc(s/pi) - 2*cos(s**2/(2*pi))*fresnels(s/pi)
         + sqrt(2)*cos(s**2/(2*pi) + pi/4))/(2*s), 0, True))
 
+    # Matrix tests
     Mt = Matrix([[exp(t), t*exp(-t)], [t*exp(-t), exp(t)]])
     Ms = Matrix([[    1/(s - 1), (s + 1)**(-2)],
                  [(s + 1)**(-2),     1/(s - 1)]])
@@ -560,35 +759,18 @@ def test_laplace_transform():
     # The default behaviour for Laplace tranform of a Matrix returns a Matrix
     # of Tuples and is deprecated:
     with warns_deprecated_sympy():
-        Ms_conds = Matrix([[(1/(s - 1), 0, Abs(s) > 1), ((s + 1)**(-2),
-            0, True)], [((s + 1)**(-2), 0, True), (1/(s - 1), 0, Abs(s) > 1)]])
+        Ms_conds = Matrix([[(1/(s - 1), 1, True), ((s + 1)**(-2),
+            -1, True)], [((s + 1)**(-2), -1, True), (1/(s - 1), 1, True)]])
     with warns_deprecated_sympy():
         assert LT(Mt, t, s) == Ms_conds
-
     # The new behavior is to return a tuple of a Matrix and the convergence
     # conditions for the matrix as a whole:
-    assert LT(Mt, t, s, legacy_matrix=False) == (Ms, 0, Abs(s) > 1)
-
+    assert LT(Mt, t, s, legacy_matrix=False) == (Ms, 1, True)
     # With noconds=True the transformed matrix is returned without conditions
     # either way:
     assert LT(Mt, t, s, noconds=True) == Ms
     assert LT(Mt, t, s, legacy_matrix=False, noconds=True) == Ms
 
-
-@slow
-def test_issue_8368t_7173():
-    LT = laplace_transform
-    # hyperbolic
-    assert LT(sinh(x), x, s) == (1/(s**2 - 1), 0, abs(s) > 1)
-    assert LT(cosh(x), x, s) == (s/(s**2 - 1), -oo, s**2 > 1)
-    assert LT(sinh(x + 3), x, s) == (
-        (-s + (s + 1)*exp(6) + 1)*exp(-3)/(s - 1)/(s + 1)/2, 0, Abs(s) > 1)
-    assert LT(sinh(x)*cosh(x), x, s) == (
-        1/(s**2 - 4), 0, Abs(s) > 2)
-    # trig (make sure they are not being rewritten in terms of exp)
-    assert LT(cos(x + 3), x, s) == ((s*cos(3) - sin(3))/(s**2 + 1), 0, True)
-
-
 @slow
 def test_inverse_laplace_transform():
     from sympy.core.exprtools import factor_terms
@@ -868,22 +1050,6 @@ def test_issue_8882():
         **{'as_meijerg': True, 'needeval': True}))
 
 
-def test_issue_7173():
-    from sympy.simplify.cse_main import cse
-    x0, x1, x2, x3 = symbols('x:4')
-    ans = laplace_transform(sinh(a*x)*cosh(a*x), x, s)
-    r, e = cse(ans)
-    assert r == [
-        (x0, arg(a)),
-        (x1, Abs(x0)),
-        (x2, pi/2),
-        (x3, Abs(x0 + pi))]
-    assert e == [
-        a/(-4*a**2 + s**2),
-        0,
-        ((x1 <= x2) | (x1 < x2)) & ((x3 <= x2) | (x3 < x2))]
-
-
 def test_issue_8514():
     from sympy.simplify.simplify import simplify
     a, b, c, = symbols('a b c', positive=True)
diff --git a/sympy/utilities/tests/test_wester.py b/sympy/utilities/tests/test_wester.py
index ea346816fb..f5c03576e6 100644
--- a/sympy/utilities/tests/test_wester.py
+++ b/sympy/utilities/tests/test_wester.py
@@ -2933,7 +2933,6 @@ def test_Y4():
     assert F == (1 - exp(-6*sqrt(s)))/s
 
 
-@XFAIL
 def test_Y5_Y6():
 # Solve y'' + y = 4 [H(t - 1) - H(t - 2)], y(0) = 1, y'(0) = 0 where H is the
 # Heaviside (unit step) function (the RHS describes a pulse of magnitude 4 and
@@ -2949,10 +2948,9 @@ def test_Y5_Y6():
                                 + y(t)
                                 - 4*(Heaviside(t - 1)
                                 - Heaviside(t - 2)), t, s)
-    # Laplace transform for diff() not calculated
-    # https://github.com/sympy/sympy/issues/7176
-    assert (F == s**2*LaplaceTransform(y(t), t, s) - s
-            + LaplaceTransform(y(t), t, s) - 4*exp(-s)/s + 4*exp(-2*s)/s)
+    assert (F == s**2*LaplaceTransform(y(t), t, s) - s*y(0) +
+            LaplaceTransform(y(t), t, s) - Subs(Derivative(y(t), t), t, 0) -
+            4*exp(-s)/s + 4*exp(-2*s)/s)
 # TODO implement second part of test case
 # Now, solve for Y(s) and then take the inverse Laplace transform
 #   => Y(s) = s/(s^2 + 1) + 4 [1/s - s/(s^2 + 1)] [e^(-s) - e^(-2 s)]